# Summary: Creates Appendix D

######################################################
#-------------------Appendix D ----------------------#
######################################################

options(scipen=999)
rm(list = ls())

# Load libraries
library(rddensity)
library(rdrobust)
library(lmtest)
library(sandwich)

# Set WD
# setwd('~/Dataverse/')

# Source functions
source('./Code/fun_analysis.R')

# Load Data
load(file = "./Data/d.RData")

# Period
d$period <- 'weak'
d$period[d$ANO_ELEICAO > 2007] <- 'strong'
d$period[d$ANO_ELEICAO > 2015] <- 'moderate'

# Centralized and Decentralized parties
d$centralized <- as.numeric(d$SIGLA_PARTIDO %in% c("PC do B", "PT", "PSDB", "PPS", "PFL", "DEM", "PL", "PR", "PDT"))

# Get codeyear for where a centralized and decentralized party won
centParties <- d$codeyear[(d$SIGLA_PARTIDO %in% c("PC do B", "PT", "PSDB", "PPS", "PFL", "DEM", "PL", "PR", "PDT")) & d$eleito == 1]
decentParties <- d$codeyear[!(d$SIGLA_PARTIDO %in% c("PC do B", "PT", "PSDB", "PPS", "PFL", "DEM", "PL", "PR", "PDT")) & d$eleito == 1]

# Keep where centralized parties won
d_d <- subset(d, codeyear %in% decentParties)
d <- subset(d, codeyear %in% centParties)

# Keep only observations from municipalities where both parties had male and female candidates
d1 <- subset(d, nonZeroDF == 2)
d_d1 <- subset(d_d, nonZeroDF == 2)
length(unique(d1$codeyear))


###########################################################################
#-------------------Appendix D.1: Descriptive Stats ----------------------#
###########################################################################

# TABLE D.8
sink("./AppendixDResults/tableD8.txt")
print("TABLE D.8: Descriptive Statistics - Only Centralized Parties")
summary(d1$avg_male_DF[d1$period == 'weak']); length(d1$avg_male_DF[d1$period == 'weak']); sd(d1$avg_male_DF[d1$period == 'weak'])
summary(d1$avg_fem_DF[d1$period == 'weak']); length(d1$avg_fem_DF[d1$period == 'weak']); sd(d1$avg_fem_DF[d1$period == 'weak'])
summary(d1$avg_male_DF[d1$period == 'strong']); length(d1$avg_male_DF[d1$period == 'strong']); sd(d1$avg_male_DF[d1$period == 'strong'])
summary(d1$avg_fem_DF[d1$period == 'strong']); length(d1$avg_fem_DF[d1$period == 'strong']); sd(d1$avg_fem_DF[d1$period == 'strong'])
summary(d1$avg_male_DF[d1$period == 'moderate']); length(d1$avg_male_DF[d1$period == 'moderate']); sd(d1$avg_male_DF[d1$period == 'moderate'])
summary(d1$avg_fem_DF[d1$period == 'moderate']); length(d1$avg_fem_DF[d1$period == 'moderate']); sd(d1$avg_fem_DF[d1$period == 'moderate'])
sink()

######################################################################
#-------------------Appendix D.2: Manipulation ----------------------#
######################################################################

run_varT1 <- d1$run_var[d1$period == 'weak']
run_varT2 <- d1$run_var[d1$period == 'strong']
run_varT3 <- d1$run_var[d1$period == 'moderate']
summary(t1 <- rddensity(X = run_varT1, c = 0)) # With massPoints
summary(t2 <- rddensity(X = run_varT2, c = 0)) # With massPoints
summary(t3 <- rddensity(X = run_varT3, c = 0)) # With massPoints
# Result: No sorting

# FIGURE D.4
# Density Plot 
pdf('./AppendixDResults/figD4.pdf', width = 8, height = 8)
par(mar = c(5, 5, 2, 2))
hist(run_varT1[abs(run_varT1)<=(13.039 - 10.873)/2 + 10.873], breaks = 50, # at the bw
	main = "", xlab = "Difference in Vote Share", cex.lab = 1.5, cex.axis = 1.5)
abline(v = 0, col = "red", lwd = 1.5)
text(x = 8.5, y = 80, labels = paste("p-value = ", round(t1$test$p_jk, 4)), cex = 1.25)
par(mar = c(5, 5, 2, 2))
hist(run_varT2[abs(run_varT2)<=(13.039 - 10.873)/2 + 10.873], breaks = 50, # at the bw
	main = "", xlab = "Difference in Vote Share", cex.lab = 1.5, cex.axis = 1.5)
abline(v = 0, col = "red", lwd = 1.5)
text(x = 8.5, y = 85, labels = paste("p-value = ", signif(t2$test$p_jk, 4)), cex = 1.25)
par(mar = c(5, 5, 2, 2))
hist(run_varT3[abs(run_varT3)<=(13.039 - 10.873)/2 + 10.873], breaks = 50, # at the bw
	main = "", xlab = "Difference in Vote Share", cex.lab = 1.5, cex.axis = 1.5)
abline(v = 0, col = "red", lwd = 1.5)
text(x = 8.5, y = 40, labels = paste("p-value = ", signif(t3$test$p_jk, 4)), cex = 1.25)
dev.off()

######################################################################
#-------------------Appendix D.3: Balance Test ----------------------#
######################################################################

################################
########## Municipal-Party Level
################################

#####--------------- Lag DV ---------------##### 
# % of vote share male
d$temp <- ifelse(is.na(d$perc_party_mun_fem_placebo1DF), 0, d$perc_party_mun_fem_placebo1DF)
d$perc_party_mun_male_placebo1DF <- d$perc_party_mun_placebo1DF - d$temp
# N of male candidates
d$temp <- ifelse(is.na(d$n_party_Femcandidates_placebo1DF), 0, d$n_party_Femcandidates_placebo1DF)
d$n_party_Malecandidates_placebo1DF <- d$n_party_candidates_placebo1DF - d$n_party_Femcandidates_placebo1DF
# Avg vote share male
d$avg_placebo1DF <- d$perc_party_mun_placebo1DF/(d$n_party_candidates_placebo1DF)
# Avg vote share male
d$avg_male_placebo1DF <- d$perc_party_mun_male_placebo1DF/(d$n_party_Malecandidates_placebo1DF)
# Avg vote share female
d$avg_fem_placebo1DF <- d$perc_party_mun_fem_placebo1DF/(d$n_party_Femcandidates_placebo1DF)

# Find municipalities that should be included (both parties had male and female candidates)
valid <- aggregate(d$n_party_Malecandidates_placebo1DF > 0 & d$n_party_Femcandidates_placebo1DF>0 & (!is.na(d$avg_male_placebo1DF) & !is.na(d$avg_fem_placebo1DF)), list(d$codeyear), sum)
valid_place <- na.omit(valid$Group.1[valid$x == 2])
dBalance <- subset(d, codeyear %in% valid_place & nonZeroDF == 2)


# Week
b1 = rdrobust(y = subset(dBalance,  period == 'weak')$avg_male_placebo1DF
		, x = subset(dBalance, period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(dBalance,  period == 'weak')$codeyear
		, all = TRUE)
summary(b1)
b2 = rdrobust(y = subset(dBalance, period == 'weak')$avg_fem_placebo1DF
		, x = subset(dBalance, period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(dBalance, period == 'weak')$codeyear
		, all = TRUE)
summary(b2)

# Strong
b3 = rdrobust(y = subset(dBalance, period == 'strong')$avg_male_placebo1DF
		, x = subset(dBalance, period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(dBalance, period == 'strong')$codeyear
		, all = TRUE)
summary(b3)
b4 = rdrobust(y = subset(dBalance,  period == 'strong')$avg_fem_placebo1DF
		, x = subset(dBalance, period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(dBalance,  period == 'strong')$codeyear
		, all = TRUE)
summary(b4)

# Moderate
b5 = rdrobust(y = subset(dBalance,  period == 'moderate')$avg_male_placebo1DF
		, x = subset(dBalance, period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(dBalance,  period == 'moderate')$codeyear
		, all = TRUE)
summary(b5)
b6 = rdrobust(y = subset(dBalance,  period == 'moderate')$avg_fem_placebo1DF
		, x = subset(dBalance,period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(dBalance, period == 'moderate')$codeyear
		, all = TRUE)
summary(b6)


#####--------------- Lag DV State Deputy---------------##### 
# % of vote share male
d$temp <- ifelse(is.na(d$perc_party_mun_fem_placebo1DE), 0, d$perc_party_mun_fem_placebo1DE)
d$perc_party_mun_male_placebo1DE <- d$perc_party_mun_placebo1DE - d$temp
# N of male candidates
d$temp <- ifelse(is.na(d$n_party_Femcandidates_placebo1DE), 0, d$n_party_Femcandidates_placebo1DE)
d$n_party_Malecandidates_placebo1DE <- d$n_party_candidates_placebo1DE - d$n_party_Femcandidates_placebo1DE
# Avg vote share male
d$avg_placebo1DE <- d$perc_party_mun_placebo1DE/(d$n_party_candidates_placebo1DE)
# Avg vote share male
d$avg_male_placebo1DE <- d$perc_party_mun_male_placebo1DE/(d$n_party_Malecandidates_placebo1DE)
# Avg vote share female
d$avg_fem_placebo1DE <- d$perc_party_mun_fem_placebo1DE/(d$n_party_Femcandidates_placebo1DE)

# Find municipalities that should be included (both parties had male and female candidates)
valid <- aggregate(d$n_party_Malecandidates_placebo1DE > 0 & d$n_party_Femcandidates_placebo1DE>0 & (!is.na(d$avg_male_placebo1DE) & !is.na(d$avg_fem_placebo1DE)), list(d$codeyear), sum)
valid_place <- na.omit(valid$Group.1[valid$x == 2])
dBalance <- subset(d, codeyear %in% valid_place & nonZeroDF == 2)

# Week
b7 = rdrobust(y = subset(dBalance,  period == 'weak')$avg_male_placebo1DE
                 , x = subset(dBalance, period == 'weak')$run_var
                 , p = 1
                 , c = 0	
                 , cluster = subset(dBalance,  period == 'weak')$codeyear
                 , all = TRUE)
summary(b7)
b8 = rdrobust(y = subset(dBalance, period == 'weak')$avg_fem_placebo1DE
                 , x = subset(dBalance, period == 'weak')$run_var
                 , p = 1
                 , c = 0	
                 , cluster = subset(dBalance, period == 'weak')$codeyear
                 , all = TRUE)
summary(b8)

# Strong
b9 = rdrobust(y = subset(dBalance, period == 'strong')$avg_male_placebo1DE
                 , x = subset(dBalance, period == 'strong')$run_var
                 , p = 1
                 , c = 0	
                 , cluster = subset(dBalance, period == 'strong')$codeyear
                 , all = TRUE)
summary(b9)
b10 = rdrobust(y = subset(dBalance,  period == 'strong')$avg_fem_placebo1DE
                 , x = subset(dBalance, period == 'strong')$run_var
                 , p = 1
                 , c = 0	
                 , cluster = subset(dBalance,  period == 'strong')$codeyear
                 , all = TRUE)
summary(b10)

# Moderate
b11 = rdrobust(y = subset(dBalance,  period == 'moderate')$avg_male_placebo1DE
                 , x = subset(dBalance, period == 'moderate')$run_var
                 , p = 1
                 , c = 0	
                 , cluster = subset(dBalance,  period == 'moderate')$codeyear
                 , all = TRUE)
summary(b11)
b12 = rdrobust(y = subset(dBalance,  period == 'moderate')$avg_fem_placebo1DE
                 , x = subset(dBalance,period == 'moderate')$run_var
                 , p = 1
                 , c = 0	
                 , cluster = subset(dBalance, period == 'moderate')$codeyear
                 , all = TRUE)
summary(b12)


#####--------------- Local Council---------------##### 

load(file = "./Data/d_balance.RData")

# Combine datasets
d <- merge(d, d_balance, by = c('codeyear', 'NUMERO_PARTIDO'), all.x = TRUE)

#### Variable for balance test
# Perct of vote share
d$perc_party_mun_placebo_Ver <- (d$total_votes_party_vereador/d$total_votes_district_vereador) * 100
d$perc_party_mun_fem_placebo_Ver <- (d$total_votes_party_fem_vereador/d$total_votes_district_vereador) * 100
d$temp <- ifelse(is.na(d$perc_party_mun_fem_placebo_Ver), 0, d$perc_party_mun_fem_placebo_Ver)
d$perc_party_mun_male_placebo_Ver <- d$perc_party_mun_placebo_Ver - d$temp 
# N of Candidates
d$temp <- ifelse(is.na(d$n_candidates_fem_vereador), 0, d$n_candidates_fem_vereador)
d$n_candidates_male_vereador <- d$n_candidates_vereador - d$temp
# Avg Vote Sahre
d$avg_vs_ver <- (d$perc_party_mun_placebo_Ver/d$n_candidates_vereador)
d$avg_vs_male_ver <- (d$perc_party_mun_male_placebo_Ver/d$n_candidates_male_vereador) 
d$avg_vs_fem_ver <- (d$perc_party_mun_fem_placebo_Ver/d$n_candidates_fem_vereador) 

# Keep municipalities that had men and women running  
valid <- aggregate((d$n_candidates_vereador - d$n_candidates_fem_vereador) > 0 & d$n_candidates_fem_vereador>0 & (!is.na(d$avg_vs_male_ver) & !is.na(d$avg_vs_fem_ver)), list(d$codeyear), sum)
valid_place <- na.omit(valid$Group.1[valid$x == 2])
dBalance <- subset(d, codeyear %in% valid_place & nonZeroDF == 2)


# Week
b13 = rdrobust(y = subset(dBalance,  period == 'weak')$avg_vs_male_ver
                 , x = subset(dBalance, period == 'weak')$run_var
                 , p = 1
                 , c = 0	
                 , cluster = subset(dBalance,  period == 'weak')$codeyear
                 , all = TRUE)
summary(b13)
b14 = rdrobust(y = subset(dBalance, period == 'weak')$avg_vs_fem_ver
                 , x = subset(dBalance, period == 'weak')$run_var
                 , p = 1
                 , c = 0	
                 , cluster = subset(dBalance, period == 'weak')$codeyear
                 , all = TRUE)
summary(b14)

# Strong
b15 = rdrobust(y = subset(dBalance, period == 'strong')$avg_vs_male_ver
                 , x = subset(dBalance, period == 'strong')$run_var
                 , p = 1
                 , c = 0	
                 , cluster = subset(dBalance, period == 'strong')$codeyear
                 , all = TRUE)
summary(b15)
b16 = rdrobust(y = subset(dBalance,  period == 'strong')$avg_vs_fem_ver
                 , x = subset(dBalance, period == 'strong')$run_var
                 , p = 1
                 , c = 0	
                 , cluster = subset(dBalance,  period == 'strong')$codeyear
                 , all = TRUE)
summary(b16)

# Moderate
b17 = rdrobust(y = subset(dBalance,  period == 'moderate')$avg_vs_male_ver
                 , x = subset(dBalance, period == 'moderate')$run_var
                 , p = 1
                 , c = 0	
                 , cluster = subset(dBalance,  period == 'moderate')$codeyear
                 , all = TRUE)
summary(b17)
b18 = rdrobust(y = subset(dBalance,  period == 'moderate')$avg_vs_fem_ver
                 , x = subset(dBalance,period == 'moderate')$run_var
                 , p = 1
                 , c = 0	
                 , cluster = subset(dBalance, period == 'moderate')$codeyear
                 , all = TRUE)
summary(b18)

# TABLE D.9
ordering = c(1,7,13,2,8,14,3,9,15,4,10,16,5,11,17,6,12,18)
eval(parse(text = paste0("Estimate = rbind(", paste0(paste0("b", ordering, "[['coef']][3]"), collapse = ","), ")") 
))
eval(parse(text = paste0("p_value = rbind(", paste0(paste0("b", ordering, "[['pv']][3]"), collapse = ","), ")") 
))
eval(parse(text = paste0("h = rbind(", paste0(paste0("b", ordering, "[['bws']][1]"), collapse = ","), ")") 
))
eval(parse(text = paste0("n_co = rbind(", paste0(paste0("b", ordering, "[['N_h']][1]"), collapse = ","), ")") 
))
eval(parse(text = paste0("n_tr = rbind(", paste0(paste0("b", ordering, "[['N_h']][2]"), collapse = ","), ")") 
))
table_d9 = cbind.data.frame(Estimate,p_value,h,n_co,n_tr)

sink("./AppendixDResults/tableD9.txt")
print("TABLE D.9: Municipal-Party Level Tests, Only Centralized Party Sample by Party Regimes")
table_d9
sink()

################################
########## Individual Level 
################################

load(file = "./Data/balance_ind.RData")

### Add variables
# Note: Unfortunately, we had to use a new version of TSE data in order to create
# the variables for the balance tests. As a result, our sample is slightly different
# Our results are robust to the use of this different sample. Also, these data are
# only available for elections after 1996.
dBalance <- merge(d, balance_ind, 
                  by.x = c('ANO_ELEICAO.x', 'CODIGO_MUNICIPIO.x', 
                           'NUMERO_PARTIDO'),
                  by.y = c('ANO_ELEICAO', 'SG_UE', 
                           'NR_PARTIDO'))
dBalance$women <- as.numeric(dBalance$CODIGO_SEXO == 4)

# Keep only municipalities that both candidates are available
dBalance$x <- 1
valid <- aggregate(dBalance$x, list(dBalance$codeyear), sum)
valid_place <- na.omit(valid$Group.1[valid$x == 2])
dBalance <- subset(dBalance, codeyear %in% valid_place & nonZeroDF == 2)


# TABLE D.10
######### Weak Period
b1 = rdrobust(y = subset(dBalance,  period == 'weak')$business_occ,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b1)
b2 = rdrobust(y = subset(dBalance,  period == 'weak')$health_occ,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b2)
b3 = rdrobust(y = subset(dBalance,  period == 'weak')$education_occ,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b3)
b4 = rdrobust(y = subset(dBalance,  period == 'weak')$merchant_occ,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b4)
b5 = rdrobust(y = subset(dBalance,  period == 'weak')$farmer_occ,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b5)
b6 = rdrobust(y = subset(dBalance,  period == 'weak')$pubservant_occ,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b6)
b7 = rdrobust(y = subset(dBalance,  period == 'weak')$laywer_pcc,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b7)
b8 = rdrobust(y = subset(dBalance,  period == 'weak')$incumbent,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b8)
b9 = rdrobust(y = subset(dBalance,  period == 'weak')$held_office,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b9)
b10 = rdrobust(y = subset(dBalance,  period == 'weak')$ran_office,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b10)

# Age is missing for one candidate in one municipality, we remove it to keep the 
# data balanced. 
keep <- names(table(dBalance$codeyear[!is.na(dBalance$age)])[table(dBalance$codeyear[!is.na(dBalance$age)])==2])
b11 = rdrobust(y = subset(dBalance,  period == 'weak' & codeyear %in% keep)$age,
                 x = subset(dBalance,  period == 'weak' & codeyear %in% keep)$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak' & codeyear %in% keep)$codeyear)
summary(b11)
b12 = rdrobust(y = subset(dBalance,  period == 'weak')$married,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b12)
b13 = rdrobust(y = subset(dBalance,  period == 'weak')$highschool,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b13)
b14 = rdrobust(y = subset(dBalance,  period == 'weak')$somecollege,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b14)
b15 = rdrobust(y = subset(dBalance,  period == 'weak')$women,
                 x = subset(dBalance,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'weak')$codeyear)
summary(b15)

######### Strong Period
b16 = rdrobust(y = subset(dBalance,  period == 'strong')$business_occ,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b16)
b17 = rdrobust(y = subset(dBalance,  period == 'strong')$health_occ,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b17)
b18 = rdrobust(y = subset(dBalance,  period == 'strong')$education_occ,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b18)
b19 = rdrobust(y = subset(dBalance,  period == 'strong')$merchant_occ,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b19)
b20 = rdrobust(y = subset(dBalance,  period == 'strong')$farmer_occ,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b20)
b21 = rdrobust(y = subset(dBalance,  period == 'strong')$pubservant_occ,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b21)
b22 = rdrobust(y = subset(dBalance,  period == 'strong')$laywer_pcc,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b22)
b23 = rdrobust(y = subset(dBalance,  period == 'strong')$incumbent,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b23)
b24 = rdrobust(y = subset(dBalance,  period == 'strong')$held_office,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b24)
b25 = rdrobust(y = subset(dBalance,  period == 'strong')$ran_office,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b25)
b26 = rdrobust(y = subset(dBalance,  period == 'strong')$age,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b26)
b27 = rdrobust(y = subset(dBalance,  period == 'strong')$married,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b27)
b28 = rdrobust(y = subset(dBalance,  period == 'strong')$highschool,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b28)
b29 = rdrobust(y = subset(dBalance,  period == 'strong')$somecollege,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b29)
b30 = rdrobust(y = subset(dBalance,  period == 'strong')$women,
                 x = subset(dBalance,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'strong')$codeyear)
summary(b30)

######### Moderate Period
b31 = rdrobust(y = subset(dBalance,  period == 'moderate')$business_occ,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b31)
b32 = rdrobust(y = subset(dBalance,  period == 'moderate')$health_occ,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b32)
b33 = rdrobust(y = subset(dBalance,  period == 'moderate')$education_occ,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b33)
b34 = rdrobust(y = subset(dBalance,  period == 'moderate')$merchant_occ,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b34)
b35 = rdrobust(y = subset(dBalance,  period == 'moderate')$farmer_occ,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b35)
b36 = rdrobust(y = subset(dBalance,  period == 'moderate')$pubservant_occ,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b36)
b37 = rdrobust(y = subset(dBalance,  period == 'moderate')$laywer_pcc,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b37)
b38 = rdrobust(y = subset(dBalance,  period == 'moderate')$incumbent,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b38)
b39 = rdrobust(y = subset(dBalance,  period == 'moderate')$held_office,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b39)
b40 = rdrobust(y = subset(dBalance,  period == 'moderate')$ran_office,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b40)
b41 = rdrobust(y = subset(dBalance,  period == 'moderate')$age,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b41)
b42 = rdrobust(y = subset(dBalance,  period == 'moderate')$married,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b42)
b43 = rdrobust(y = subset(dBalance,  period == 'moderate')$highschool,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b43)
b44 = rdrobust(y = subset(dBalance,  period == 'moderate')$somecollege,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b44)
b45 = rdrobust(y = subset(dBalance,  period == 'moderate')$women,
                 x = subset(dBalance,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance,  period == 'moderate')$codeyear)
summary(b45)

# Campaign Expeditures are only available after 2000 and many candidates do not report
# So, we will do a separate analysis
load(file = './Data/DespesaMunicipal.RData')

# Add Data
dBalance2 <- merge(dBalance, despesa, 
                   by.x = c('CODIGO_MUNICIPIO.x', 'ANO_ELEICAO.x', 'NUMERO_PARTIDO'), 
                   by.y = c('SG_UE', 'ANO_ELEICAO', 'NR_PART'))
# Keep only municipalities that both candidates are available
dBalance2$x <- 1
valid <- aggregate(dBalance2$x, list(dBalance2$codeyear), sum)
valid_place <- na.omit(valid$Group.1[valid$x == 2])
dBalance2 <- subset(dBalance2, codeyear %in% valid_place & nonZeroDF == 2)

# Rescale by 1000
dBalance2$DESPESA <- dBalance2$DESPESA/1000

b46 = rdrobust(y = subset(dBalance2,  period == 'weak')$DESPESA,
                 x = subset(dBalance2,  period == 'weak')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance2,  period == 'weak')$codeyear)
summary(b46)
b47 = rdrobust(y = subset(dBalance2,  period == 'strong')$DESPESA,
                 x = subset(dBalance2,  period == 'strong')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance2,  period == 'strong')$codeyear)
summary(b47)
b48 = rdrobust(y = subset(dBalance2,  period == 'moderate')$DESPESA,
                 x = subset(dBalance2,  period == 'moderate')$run_var,
                 p = 1,
                 c = 0,
                 all = TRUE,
                 cluster = subset(dBalance2,  period == 'moderate')$codeyear)
summary(b48)

# TABLE D.10
ordering = c(1:48)
ordering = c(1:15,46,16:30,47,31:45,48)

eval(parse(text = paste0("Estimate = rbind(", paste0(paste0("b", ordering, "[['coef']][3]"), collapse = ","), ")") 
))
eval(parse(text = paste0("p_value = rbind(", paste0(paste0("b", ordering, "[['pv']][3]"), collapse = ","), ")") 
))
eval(parse(text = paste0("h = rbind(", paste0(paste0("b", ordering, "[['bws']][1]"), collapse = ","), ")") 
))
eval(parse(text = paste0("n_co = rbind(", paste0(paste0("b", ordering, "[['N_h']][1]"), collapse = ","), ")") 
))
eval(parse(text = paste0("n_tr = rbind(", paste0(paste0("b", ordering, "[['N_h']][2]"), collapse = ","), ")") 
))
table_d10 = cbind.data.frame(Estimate,p_value,h,n_co,n_tr)

sink("./AppendixDResults/tableD10.txt")
print("TABLE D.10: Candidate Level Tests, Only Centralized Party Sample by Party Regimes")
table_d10
sink()


### Run Models with Controls
# TABLE D.11
# Weak
out.cent.men.opt.weak <- rdrobust(y = subset(dBalance, period == 'weak')$avg_male_DF
                                  , x = subset(dBalance, period == 'weak')$run_var
                                  , covs = cbind(subset(dBalance, period == 'weak')$health_occ, 
                                                 subset(dBalance, period == 'weak')$married)
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(dBalance, period == 'weak')$codeyear
                                  , all = TRUE)
out.cent.women.opt.weak <- rdrobust(y = subset(dBalance, period == 'weak')$avg_fem_DF
                                    , x = subset(dBalance, period == 'weak')$run_var
                                    , covs = cbind(subset(dBalance, period == 'weak')$health_occ, 
                                                   subset(dBalance, period == 'weak')$married)                                 
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(dBalance, period == 'weak')$codeyear
                                    , all = TRUE)
summary(out.cent.men.opt.weak)
summary(out.cent.women.opt.weak)

# Strong
# Find municipalities that should be included (both parties had male and female candidates in terms of state deputies)
valid <- aggregate(dBalance$n_party_Malecandidates_placebo1DE > 0 & dBalance$n_party_Femcandidates_placebo1DE>0 & 
                     (!is.na(dBalance$avg_male_placebo1DE) & !is.na(dBalance$avg_fem_placebo1DE)) &
                     !is.na(subset(dBalance)$health_occ) &
                     !is.na(subset(dBalance)$age) 
                   , list(dBalance$codeyear), sum)
valid_place <- na.omit(valid$Group.1[valid$x == 2])
dBalance2 <- subset(dBalance, codeyear %in% valid_place & nonZeroDF == 2)

out.cent.men.opt.strong <- rdrobust(y = subset(dBalance2, period == 'strong')$avg_male_DF
                                    , x = subset(dBalance2, period == 'strong')$run_var
                                    , covs = cbind(subset(dBalance2, period == 'strong')$age,
                                                   subset(dBalance2, period == 'strong')$avg_male_placebo1DE)                                    
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(dBalance2, period == 'strong')$codeyear
                                    , all = TRUE)
out.cent.women.opt.strong <- rdrobust(y = subset(dBalance2, period == 'strong')$avg_fem_DF
                                      , x = subset(dBalance2, period == 'strong')$run_var
                                      , covs = cbind(subset(dBalance2, period == 'strong')$age,
                                                     subset(dBalance2, period == 'strong')$avg_male_placebo1DE)                                        
                                      , p = 1
                                      , c = 0	
                                      , cluster = subset(dBalance2, period == 'strong')$codeyear
                                      , all = TRUE)
summary(out.cent.men.opt.strong)
summary(out.cent.women.opt.strong)

# Strong
out.cent.men.opt.strong2 <- rdrobust(y = subset(dBalance, period == 'strong')$avg_male_DF
                                      , x = subset(dBalance, period == 'strong')$run_var
                                      , covs = cbind(subset(dBalance, period == 'strong')$age)                                         
                                      , p = 1
                                      , c = 0	
                                      , cluster = subset(dBalance, period == 'strong')$codeyear
                                      , all = TRUE)
out.cent.women.opt.strong2 <- rdrobust(y = subset(dBalance, period == 'strong')$avg_fem_DF
                                        , x = subset(dBalance, period == 'strong')$run_var
                                        , covs = cbind(subset(dBalance, period == 'strong')$age)                                        
                                        , p = 1
                                        , c = 0	
                                        , cluster = subset(dBalance, period == 'strong')$codeyear
                                        , all = TRUE)
summary(out.cent.men.opt.strong2)
summary(out.cent.women.opt.strong2)


# Moderate
out.cent.men.opt.moderate <- rdrobust(y = subset(dBalance, period == 'moderate')$avg_male_DF
                                      , x = subset(dBalance, period == 'moderate')$run_var
                                      , covs = cbind(subset(dBalance, period == 'moderate')$education_occ, 
                                                     subset(dBalance, period == 'moderate')$somecollege)                                         
                                      , p = 1
                                      , c = 0	
                                      , cluster = subset(dBalance, period == 'moderate')$codeyear
                                      , all = TRUE)
out.cent.women.opt.moderate <- rdrobust(y = subset(dBalance, period == 'moderate')$avg_fem_DF
                                        , x = subset(dBalance, period == 'moderate')$run_var
                                        , covs = cbind(subset(dBalance, period == 'moderate')$education_occ, 
                                                       subset(dBalance, period == 'moderate')$somecollege)                                        
                                        , p = 1
                                        , c = 0	
                                        , cluster = subset(dBalance, period == 'moderate')$codeyear
                                        , all = TRUE)
summary(out.cent.men.opt.moderate)
summary(out.cent.women.opt.moderate)

# Calculate Differences 
zOutWeak <- zDif(beta1 = out.cent.men.opt.weak$coef[3,], beta2 = out.cent.women.opt.weak$coef[3,]
                 , se1 = out.cent.men.opt.weak$se[3,], se2 = out.cent.women.opt.weak$se[3,])
zOutStrong <- zDif(beta1 = out.cent.men.opt.strong$coef[3,], beta2 = out.cent.women.opt.strong$coef[3,]
                   , se1 = out.cent.men.opt.strong$se[3,], se2 = out.cent.women.opt.strong$se[3,])
zOutStrong2 <- zDif(beta1 = out.cent.men.opt.strong2$coef[3,], beta2 = out.cent.women.opt.strong2$coef[3,]
                   , se1 = out.cent.men.opt.strong2$se[3,], se2 = out.cent.women.opt.strong2$se[3,])
zOutModerate <- zDif(beta1 = out.cent.men.opt.moderate$coef[3,], beta2 = out.cent.women.opt.moderate$coef[3,]
                     , se1 = out.cent.men.opt.moderate$se[3,], se2 = out.cent.women.opt.moderate$se[3,])

zOutWeak
(pnorm(zOutWeak$zscore)) * 2 

zOutStrong
(1-pnorm(zOutStrong$zscore)) * 2 

zOutStrong2
(1-pnorm(zOutStrong2$zscore)) * 2 

zOutModerate
(1-pnorm(zOutModerate$zscore)) * 2 

# TABLE D.11
Estimate = rbind(out.cent.men.opt.weak$coef[3],out.cent.women.opt.weak$coef[3],zOutWeak[1],
                 out.cent.men.opt.strong$coef[3],out.cent.women.opt.strong$coef[3],zOutStrong[1],
                 out.cent.men.opt.strong2$coef[3],out.cent.women.opt.strong2$coef[3],zOutStrong2[1],
                 out.cent.men.opt.moderate$coef[3],out.cent.women.opt.moderate$coef[3],zOutModerate[1])
p_value = rbind(out.cent.men.opt.weak$pv[3],out.cent.women.opt.weak$pv[3],(pnorm(zOutWeak$zscore)) * 2,
                 out.cent.men.opt.strong$pv[3],out.cent.women.opt.strong$pv[3],(1-pnorm(zOutStrong$zscore)) * 2 ,
                 out.cent.men.opt.strong2$pv[3],out.cent.women.opt.strong2$pv[3],(1-pnorm(zOutStrong2$zscore)) * 2 ,
                 out.cent.men.opt.moderate$pv[3],out.cent.women.opt.moderate$pv[3],(1-pnorm(zOutModerate$zscore)) * 2 )
h = rbind(out.cent.men.opt.weak$bws[1],out.cent.women.opt.weak$bws[1],NA,
          out.cent.men.opt.strong$bws[1],out.cent.women.opt.strong$bws[1],NA,
          out.cent.men.opt.strong2$bws[1],out.cent.women.opt.strong2$bws[1],NA,
          out.cent.men.opt.moderate$bws[1],out.cent.women.opt.moderate$bws[1],NA)
n_co = rbind(out.cent.men.opt.weak$N_h[1],out.cent.women.opt.weak$N_h[1],NA,
             out.cent.men.opt.strong$N_h[1],out.cent.women.opt.strong$N_h[1],NA,
             out.cent.men.opt.strong2$N_h[1],out.cent.women.opt.strong2$N_h[1],NA,
             out.cent.men.opt.moderate$N_h[1],out.cent.women.opt.moderate$N_h[1],NA)
n_tr = rbind(out.cent.men.opt.weak$N_h[2],out.cent.women.opt.weak$N_h[2],NA,
             out.cent.men.opt.strong$N_h[2],out.cent.women.opt.strong$N_h[2],NA,
             out.cent.men.opt.strong2$N_h[2],out.cent.women.opt.strong2$N_h[2],NA,
             out.cent.men.opt.moderate$N_h[2],out.cent.women.opt.moderate$N_h[2],NA)
table_d11 = cbind.data.frame(Estimate,p_value,h,n_co,n_tr)

sink("./AppendixDResults/tableD11.txt")
print("TABLE D.11: RD Effect of Co-Partisan Mayor on Men and Women Candidates Vote Share—
Brazilian Federal Elections, 2000-2018 (Including Covariates that are imbalanced)")
table_d11
sink()

######################################################################
#-------------------Appendix D.4: Sensitivity ----------------------#
######################################################################

# FIGURE D.5
genSensitivityPeriodH1(yearElection = c(1996, 2000, 2004), yearPlot = 'Weak', file = "./AppendixDResults/figD5.pdf")

# FIGURE D.6
genSensitivityPeriodH1(yearElection =  c(2008, 2012), yearPlot = 'Strong', file = "./AppendixDResults/figD6.pdf")

# FIGURE D.7
genSensitivityPeriodH1(yearElection = 2016, yearPlot = 'Moderate', file = "./AppendixDResults/figD7.pdf")


###############################################################################
#-------------------Appendix D.5: Different Polynomials ----------------------#
###############################################################################

# TABLE D.12
sink("./AppendixDResults/tableD12.txt")
print("TABLE D.12: RD Effect of Co-Partisan Mayor on Average Vote Share (Different Polynomials)—
Brazilian Federal Elections, 1998-2016")

polyH1(yearElection = c(1996, 2000, 2004), p = 2)
polyH1(yearElection = c(1996, 2000, 2004), p = 3)
polyH1(yearElection = c(1996, 2000, 2004), p = 4)

polyH1(yearElection = c(2008, 2012), p = 2)
polyH1(yearElection = c(2008, 2012), p = 3)
polyH1(yearElection = c(2008, 2012), p = 4)

polyH1(yearElection = 2016, p = 2)
polyH1(yearElection = 2016, p = 3)
polyH1(yearElection = 2016, p = 4)

sink()


########################################################################
#-------------------Appendix D.6: Previously Ran ----------------------#
########################################################################

# N of male candidates
d1$n_party_candidates_male_candBefore_DF <- d1$n_party_candidates_candBefore_DF - d1$n_party_Femcandidates_candBefore_DF
valid <- aggregate(d1$n_party_candidates_male_candBefore_DF > 0 & d1$n_party_Femcandidates_candBefore_DF>0, list(d1$codeyear), sum)
valid_ranB <- na.omit(valid$Group.1[valid$x == 2])
# % of vote share male
d1$perc_party_mun_male_candBeforeDF <- d1$perc_party_mun_candBefore_DF - d1$perc_party_mun_fem_candBefore_DF
# Avg vote share all
d1$avg_candBefore_DF <- d1$perc_party_mun_candBefore_DF/d1$n_party_candidates_candBefore_DF
# Avg vote share male
d1$avg_male_candBefore_DF <- d1$perc_party_mun_male_candBeforeDF/(d1$n_party_candidates_candBefore_DF - d1$n_party_Femcandidates_candBefore_DF)
# Avg vote share female
d1$avg_fem_candBefore_DF <- d1$perc_party_mun_fem_candBefore_DF/(d1$n_party_Femcandidates_candBefore_DF)
# Subset data
d2 <- subset(d1, codeyear %in% valid_ranB)

# Weak
out.cent.men.opt.weak <- rdrobust(y = subset(d2, codeyear %in% centParties  & period == 'weak')$avg_male_candBefore_DF
		, x = subset(d2, codeyear %in% centParties  & period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties  & period == 'weak')$codeyear
		, all = TRUE)
out.cent.women.opt.weak <- rdrobust(y = subset(d2, codeyear %in% centParties  & period == 'weak')$avg_fem_candBefore_DF
		, x = subset(d2, codeyear %in% centParties  & period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties  & period == 'weak')$codeyear
		, all = TRUE)

# Strong
out.cent.men.opt.strong <- rdrobust(y = subset(d2, codeyear %in% centParties  & period == 'strong')$avg_male_candBefore_DF
		, x = subset(d2, codeyear %in% centParties  & period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties  & period == 'strong')$codeyear
		, all = TRUE)
out.cent.women.opt.strong <- rdrobust(y = subset(d2, codeyear %in% centParties  & period == 'strong')$avg_fem_candBefore_DF
		, x = subset(d2, codeyear %in% centParties  & period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties  & period == 'strong')$codeyear
		, all = TRUE)

# Moderate
out.cent.men.opt.moderate <- rdrobust(y = subset(d2, codeyear %in% centParties  & period == 'moderate')$avg_male_candBefore_DF
		, x = subset(d2, codeyear %in% centParties  & period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties  & period == 'moderate')$codeyear
		, all = TRUE)
out.cent.women.opt.moderate <- rdrobust(y = subset(d2, codeyear %in% centParties  & period == 'moderate')$avg_fem_candBefore_DF
		, x = subset(d2, codeyear %in% centParties  & period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties  & period == 'moderate')$codeyear
		, all = TRUE)

# Get point estimates, se, and ci for men
betas_men <- c(out.cent.men.opt.weak$coef[3,], out.cent.men.opt.strong$coef[3,]
	, out.cent.men.opt.moderate$coef[3,])
se_men <- c(out.cent.men.opt.weak$se[3,], out.cent.men.opt.strong$se[3,]
	, out.cent.men.opt.moderate$se[3,])
lb_men <- c(out.cent.men.opt.weak$ci[3,1], out.cent.men.opt.strong$ci[3,1]
	, out.cent.men.opt.moderate$ci[3,1])
ub_men <- c(out.cent.men.opt.weak$ci[3,2], out.cent.men.opt.strong$ci[3,2]
	, out.cent.men.opt.moderate$ci[3,2])

# Get point estimates, se, and ci for women
betas_women <- c(out.cent.women.opt.weak$coef[3,], out.cent.women.opt.strong$coef[3,]
	, out.cent.women.opt.moderate$coef[3,])
se_women <- c(out.cent.women.opt.weak$se[3,], out.cent.women.opt.strong$se[3,]
	, out.cent.women.opt.moderate$se[3,])
lb_women <- c(out.cent.women.opt.weak$ci[3,1], out.cent.women.opt.strong$ci[3,1]
	, out.cent.women.opt.moderate$ci[3,1])
ub_women <- c(out.cent.women.opt.weak$ci[3,2], out.cent.women.opt.strong$ci[3,2]
	, out.cent.women.opt.moderate$ci[3,2])

# Calculate Differences 
zOutWeak <- zDif(beta1 = out.cent.men.opt.weak$coef[3,], beta2 = out.cent.women.opt.weak$coef[3,]
	, se1 = out.cent.men.opt.weak$se[3,], se2 = out.cent.women.opt.weak$se[3,])
zOutStrong <- zDif(beta1 = out.cent.men.opt.strong$coef[3,], beta2 = out.cent.women.opt.strong$coef[3,]
	, se1 = out.cent.men.opt.strong$se[3,], se2 = out.cent.women.opt.strong$se[3,])
zOutModerate <- zDif(beta1 = out.cent.men.opt.moderate$coef[3,], beta2 = out.cent.women.opt.moderate$coef[3,]
	, se1 = out.cent.men.opt.moderate$se[3,], se2 = out.cent.women.opt.moderate$se[3,])
zOut <- list(zOutWeak, zOutStrong, zOutModerate)

# FIGURE D.8
pdf('./AppendixDResults/figD8.pdf', width = 14)
byPeriodH1(betas_men = betas_men, betas_women = betas_women
	, lb_men = lb_men, ub_men = ub_men, lb_women = lb_women
	, ub_women = ub_women, zout = zOut)
dev.off()


########################################################################
#-------------------Appendix D.7: Previously Won ----------------------#
########################################################################

# ##### Candidate who won before
d1$perc_party_mun_fem_wonBefore_DF <- (d1$total_party_mun_fem_wonBefore_DF/d1$total_mun_wonBefore_DF) * 100
# # % of vote share male
d1$perc_party_mun_male_wonBefore_DF <- d1$perc_party_mun_wonBefore_DF - d1$perc_party_mun_fem_wonBefore_DF
# # N of male candidates
d1$n_party_candidates_male_wonBefore_DF <- d1$n_party_candidates_wonBefore_DF - d1$n_party_Femcandidates_wonBefore_DF
# # Avg vote share 
d1$avg_wonBefore_DF <- d1$perc_party_mun_wonBefore_DF/d1$n_party_candidates_wonBefore_DF
# # Avg vote share male
d1$avg_male_wonBefore_DF <- d1$perc_party_mun_male_wonBefore_DF/d1$n_party_candidates_male_wonBefore_DF
# # Avg vote share female
d1$avg_fem_wonBefore_DF <- d1$perc_party_mun_fem_wonBefore_DF/(d1$n_party_Femcandidates_wonBefore_DF)

# # Find municipalities that should be included
valid <- aggregate(d1$n_party_candidates_male_wonBefore_DF > 0 & d1$n_party_Femcandidates_wonBefore_DF>0, list(d1$codeyear), sum)
valid_wonB <- na.omit(valid$Group.1[valid$x == 2])

# Subset data
d2 <- subset(d1, codeyear %in% valid_wonB)

# Weak
out.cent.men.opt.weak <- rdrobust(y = subset(d2, codeyear %in% centParties  & period == 'weak')$avg_male_wonBefore_DF
		, x = subset(d2, codeyear %in% centParties  & period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties  & period == 'weak')$codeyear
		, all = TRUE)
out.cent.women.opt.weak <- rdrobust(y = subset(d2, codeyear %in% centParties  & period == 'weak')$avg_fem_wonBefore_DF
		, x = subset(d2, codeyear %in% centParties  & period == 'weak')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties  & period == 'weak')$codeyear
		, all = TRUE)

# Strong
out.cent.men.opt.strong <- rdrobust(y = subset(d2, codeyear %in% centParties  & period == 'strong')$avg_male_wonBefore_DF
		, x = subset(d2, codeyear %in% centParties  & period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties  & period == 'strong')$codeyear
		, all = TRUE)
out.cent.women.opt.strong <- rdrobust(y = subset(d2, codeyear %in% centParties  & period == 'strong')$avg_fem_wonBefore_DF
		, x = subset(d2, codeyear %in% centParties  & period == 'strong')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties  & period == 'strong')$codeyear
		, all = TRUE)

# Moderate
out.cent.men.opt.moderate <- rdrobust(y = subset(d2, codeyear %in% centParties  & period == 'moderate')$avg_male_wonBefore_DF
		, x = subset(d2, codeyear %in% centParties  & period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties  & period == 'moderate')$codeyear
		, all = TRUE)
out.cent.women.opt.moderate <- rdrobust(y = subset(d2, codeyear %in% centParties  & period == 'moderate')$avg_fem_wonBefore_DF
		, x = subset(d2, codeyear %in% centParties  & period == 'moderate')$run_var
		, p = 1
		, c = 0	
		, cluster = subset(d2, codeyear %in% centParties  & period == 'moderate')$codeyear
		, all = TRUE)

# Get point estimates, se, and ci for men
betas_men <- c(out.cent.men.opt.weak$coef[3,], out.cent.men.opt.strong$coef[3,]
	, out.cent.men.opt.moderate$coef[3,])
se_men <- c(out.cent.men.opt.weak$se[3,], out.cent.men.opt.strong$se[3,]
	, out.cent.men.opt.moderate$se[3,])
lb_men <- c(out.cent.men.opt.weak$ci[3,1], out.cent.men.opt.strong$ci[3,1]
	, out.cent.men.opt.moderate$ci[3,1])
ub_men <- c(out.cent.men.opt.weak$ci[3,2], out.cent.men.opt.strong$ci[3,2]
	, out.cent.men.opt.moderate$ci[3,2])

# Get point estimates, se, and ci for women
betas_women <- c(out.cent.women.opt.weak$coef[3,], out.cent.women.opt.strong$coef[3,]
	, out.cent.women.opt.moderate$coef[3,])
se_women <- c(out.cent.women.opt.weak$se[3,], out.cent.women.opt.strong$se[3,]
	, out.cent.women.opt.moderate$se[3,])
lb_women <- c(out.cent.women.opt.weak$ci[3,1], out.cent.women.opt.strong$ci[3,1]
	, out.cent.women.opt.moderate$ci[3,1])
ub_women <- c(out.cent.women.opt.weak$ci[3,2], out.cent.women.opt.strong$ci[3,2]
	, out.cent.women.opt.moderate$ci[3,2])

# Calculate Differences 
zOutWeak <- zDif(beta1 = out.cent.men.opt.weak$coef[3,], beta2 = out.cent.women.opt.weak$coef[3,]
	, se1 = out.cent.men.opt.weak$se[3,], se2 = out.cent.women.opt.weak$se[3,])
zOutStrong <- zDif(beta1 = out.cent.men.opt.strong$coef[3,], beta2 = out.cent.women.opt.strong$coef[3,]
	, se1 = out.cent.men.opt.strong$se[3,], se2 = out.cent.women.opt.strong$se[3,])
zOutModerate <- zDif(beta1 = out.cent.men.opt.moderate$coef[3,], beta2 = out.cent.women.opt.moderate$coef[3,]
	, se1 = out.cent.men.opt.moderate$se[3,], se2 = out.cent.women.opt.moderate$se[3,])
zOut <- list(zOutWeak, zOutStrong, zOutModerate)


# FIGURE D.9
pdf('./AppendixDResults/figD9.pdf', width = 14)
byPeriodH1(betas_men = betas_men, betas_women = betas_women
	, lb_men = lb_men, ub_men = ub_men, lb_women = lb_women
	, ub_women = ub_women
	, zout = zOut, ylims = c(-6, 6))
dev.off()



####################################################################
#-------------------Appendix D.8: RD Linear  ----------------------#
####################################################################

# Create Dataset
dmale <- subset(d1, select = c('avg_male_DF', 'position', 'run_var', 'codeyear', 'mayorFem', 'ANO_ELEICAO'))
dmale$men <- 1
dfem <- subset(d1, select = c('avg_fem_DF', 'position', 'run_var', 'codeyear', 'mayorFem', 'ANO_ELEICAO'))
dfem$men <- 0
names(dmale)[1] <- 'avg_DF'
names(dfem)[1] <- 'avg_DF'
dall <- rbind(dmale, dfem)
dall$winner <- as.numeric(dall$run_var > 0)
d2 <- subset(dall, codeyear %in% centParties)
# Run Models
out1 <- linerEstimator(c(1996, 2000, 2004))
out2 <- linerEstimator(c(2008, 2012))
out3 <- linerEstimator(2016)

# Generate plots
# FIGURE D.10
pdf('./AppendixDResults/figD10.pdf')
plotTtest(betas = out1[[1]], lb = out1[[2]], ub = out1[[3]], dvName = 'Federal Deputy - Men', yub = 4, ylb = -4)
plotTtest(betas = out1[[4]], lb = out1[[5]], ub = out1[[6]], dvName = 'Federal Deputy - Women', yub = 4, ylb = -4)
plotTtest(betas = out1[[7]], lb = out1[[8]], ub = out1[[9]], dvName = 'Federal Deputy - Difference', yub = 4, ylb = -4)
dev.off()
# FIGURE D.11
pdf('./AppendixDResults/figD11.pdf')
plotTtest(betas = out2[[1]], lb = out2[[2]], ub = out2[[3]], dvName = 'Federal Deputy - Men', yub = 4, ylb = -4)
plotTtest(betas = out2[[4]], lb = out2[[5]], ub = out2[[6]], dvName = 'Federal Deputy - Women', yub = 4, ylb = -4)
plotTtest(betas = out2[[7]], lb = out2[[8]], ub = out2[[9]], dvName = 'Federal Deputy - Difference', yub = 4, ylb = -4)
dev.off()
# FIGURE D.12
pdf('./AppendixDResults/figD12.pdf')
plotTtest(betas = out3[[1]], lb = out3[[2]], ub = out3[[3]], dvName = 'Federal Deputy - Men', yub = 4, ylb = -4)
plotTtest(betas = out3[[4]], lb = out3[[5]], ub = out3[[6]], dvName = 'Federal Deputy - Women', yub = 4, ylb = -4)
plotTtest(betas = out3[[7]], lb = out3[[8]], ub = out3[[9]], dvName = 'Federal Deputy - Difference', yub = 4, ylb = -4)
dev.off()

########################################################################
#-------------------Appendix D.9: By Elections ----------------------#
########################################################################

# TABLE D.13
#### Descriptive Stats

sink("./AppendixDResults/tableD13.txt")
print("TABLE D.13: Descriptive Statistics - Only Centralized Parties")
summary(d1$avg_male_DF[d1$ANO_ELEICAO == 1996]); length(d1$avg_male_DF[d1$ANO_ELEICAO == 1996]); sd(d1$avg_male_DF[d1$ANO_ELEICAO == 1996])
summary(d1$avg_fem_DF[d1$ANO_ELEICAO == 1996]); length(d1$avg_fem_DF[d1$ANO_ELEICAO == 1996]); sd(d1$avg_fem_DF[d1$ANO_ELEICAO == 1996])
summary(d1$avg_male_DF[d1$ANO_ELEICAO == 2000]); length(d1$avg_male_DF[d1$ANO_ELEICAO == 2000]); sd(d1$avg_male_DF[d1$ANO_ELEICAO == 2000])
summary(d1$avg_fem_DF[d1$ANO_ELEICAO == 2000]); length(d1$avg_fem_DF[d1$ANO_ELEICAO == 2000]); sd(d1$avg_fem_DF[d1$ANO_ELEICAO == 2000])
summary(d1$avg_male_DF[d1$ANO_ELEICAO == 2004]); length(d1$avg_male_DF[d1$ANO_ELEICAO == 2004]); sd(d1$avg_male_DF[d1$ANO_ELEICAO == 2004])
summary(d1$avg_fem_DF[d1$ANO_ELEICAO == 2004]); length(d1$avg_fem_DF[d1$ANO_ELEICAO == 2004]); sd(d1$avg_fem_DF[d1$ANO_ELEICAO == 2004])
summary(d1$avg_male_DF[d1$ANO_ELEICAO == 2008]); length(d1$avg_male_DF[d1$ANO_ELEICAO == 2008]); sd(d1$avg_male_DF[d1$ANO_ELEICAO == 2008])
summary(d1$avg_fem_DF[d1$ANO_ELEICAO == 2008]); length(d1$avg_fem_DF[d1$ANO_ELEICAO == 2008]); sd(d1$avg_fem_DF[d1$ANO_ELEICAO == 2008])
summary(d1$avg_male_DF[d1$ANO_ELEICAO == 2012]); length(d1$avg_male_DF[d1$ANO_ELEICAO == 2012]); sd(d1$avg_male_DF[d1$ANO_ELEICAO == 2012])
summary(d1$avg_fem_DF[d1$ANO_ELEICAO == 2012]); length(d1$avg_fem_DF[d1$ANO_ELEICAO == 2012]); sd(d1$avg_fem_DF[d1$ANO_ELEICAO == 2012])
summary(d1$avg_male_DF[d1$ANO_ELEICAO == 2016]); length(d1$avg_male_DF[d1$ANO_ELEICAO == 2016]); sd(d1$avg_male_DF[d1$ANO_ELEICAO == 2016])
summary(d1$avg_fem_DF[d1$ANO_ELEICAO == 2016]); length(d1$avg_fem_DF[d1$ANO_ELEICAO == 2016]); sd(d1$avg_fem_DF[d1$ANO_ELEICAO == 2016])
sink()


##### Results by Election

# 1998
out.cent.men.opt.1998 <- rdrobust(y = subset(d1, ANO_ELEICAO == 1996)$avg_male_DF
                                  , x = subset(d1, codeyear %in% centParties  & ANO_ELEICAO == 1996)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d1, ANO_ELEICAO == 1996)$codeyear
                                  , all = TRUE)
out.cent.women.opt.1998 <- rdrobust(y = subset(d1, ANO_ELEICAO == 1996)$avg_fem_DF
                                    , x = subset(d1, codeyear %in% centParties  & ANO_ELEICAO == 1996)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d1, ANO_ELEICAO == 1996)$codeyear
                                    , all = TRUE)

# 2002
out.cent.men.opt.2002 <- rdrobust(y = subset(d1, ANO_ELEICAO == 2000)$avg_male_DF
                                  , x = subset(d1, ANO_ELEICAO == 2000)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d1, ANO_ELEICAO == 2000)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2002 <- rdrobust(y = subset(d1, ANO_ELEICAO == 2000)$avg_fem_DF
                                    , x = subset(d1, ANO_ELEICAO == 2000)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d1, ANO_ELEICAO == 2000)$codeyear
                                    , all = TRUE)

# 2006
out.cent.men.opt.2006 <- rdrobust(y = subset(d1, ANO_ELEICAO == 2004)$avg_male_DF
                                  , x = subset(d1, ANO_ELEICAO == 2004)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d1, ANO_ELEICAO == 2004)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2006 <- rdrobust(y = subset(d1, ANO_ELEICAO == 2004)$avg_fem_DF
                                    , x = subset(d1, ANO_ELEICAO == 2004)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d1, ANO_ELEICAO == 2004)$codeyear
                                    , all = TRUE)

# 2010
out.cent.men.opt.2010 <- rdrobust(y = subset(d1, ANO_ELEICAO == 2008)$avg_male_DF
                                  , x = subset(d1, ANO_ELEICAO == 2008)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d1, ANO_ELEICAO == 2008)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2010 <- rdrobust(y = subset(d1, ANO_ELEICAO == 2008)$avg_fem_DF
                                    , x = subset(d1, ANO_ELEICAO == 2008)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d1, ANO_ELEICAO == 2008)$codeyear
                                    , all = TRUE)

# 2014
out.cent.men.opt.2014 <- rdrobust(y = subset(d1, ANO_ELEICAO == 2012)$avg_male_DF
                                  , x = subset(d1, ANO_ELEICAO == 2012)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d1, ANO_ELEICAO == 2012)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2014 <- rdrobust(y = subset(d1, ANO_ELEICAO == 2012)$avg_fem_DF
                                    , x = subset(d1, ANO_ELEICAO == 2012)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d1, ANO_ELEICAO == 2012)$codeyear
                                    , all = TRUE)

# 2018
out.cent.men.opt.2018 <- rdrobust(y = subset(d1, ANO_ELEICAO == 2016)$avg_male_DF
                                  , x = subset(d1, ANO_ELEICAO == 2016)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d1, ANO_ELEICAO == 2016)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2018 <- rdrobust(y = subset(d1, ANO_ELEICAO == 2016)$avg_fem_DF
                                    , x = subset(d1, ANO_ELEICAO == 2016)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d1, ANO_ELEICAO == 2016)$codeyear
                                    , all = TRUE)

# Get point estimates, se, and ci for men
betas_men <- c(out.cent.men.opt.1998$coef[3,], out.cent.men.opt.2002$coef[3,]
               , out.cent.men.opt.2006$coef[3,], out.cent.men.opt.2010$coef[3,], out.cent.men.opt.2014$coef[3,]
               , out.cent.men.opt.2018$coef[3,])
se_men <- c(out.cent.men.opt.1998$se[3,], out.cent.men.opt.2002$se[3,]
            , out.cent.men.opt.2006$se[3,], out.cent.men.opt.2010$se[3,], out.cent.men.opt.2014$se[3,]
            , out.cent.men.opt.2018$se[3,])
lb_men <- c(out.cent.men.opt.1998$ci[3,1], out.cent.men.opt.2002$ci[3,1]
            , out.cent.men.opt.2006$ci[3,1], out.cent.men.opt.2010$ci[3,1], out.cent.men.opt.2014$ci[3,1]
            , out.cent.men.opt.2018$ci[3,1])
ub_men <- c(out.cent.men.opt.1998$ci[3,2], out.cent.men.opt.2002$ci[3,2]
            , out.cent.men.opt.2006$ci[3,2], out.cent.men.opt.2010$ci[3,2], out.cent.men.opt.2014$ci[3,2]
            , out.cent.men.opt.2018$ci[3,2])

# Get point estimates, se, and ci for women
betas_women <- c(out.cent.women.opt.1998$coef[3,], out.cent.women.opt.2002$coef[3,]
                 , out.cent.women.opt.2006$coef[3,], out.cent.women.opt.2010$coef[3,], out.cent.women.opt.2014$coef[3,]
                 , out.cent.women.opt.2018$coef[3,])
se_women <- c(out.cent.women.opt.1998$se[3,], out.cent.women.opt.2002$se[3,]
              , out.cent.women.opt.2006$se[3,], out.cent.women.opt.2010$se[3,], out.cent.women.opt.2014$se[3,]
              , out.cent.women.opt.2018$se[3,])
lb_women <- c(out.cent.women.opt.1998$ci[3,1], out.cent.women.opt.2002$ci[3,1]
              , out.cent.women.opt.2006$ci[3,1], out.cent.women.opt.2010$ci[3,1], out.cent.women.opt.2014$ci[3,1]
              , out.cent.women.opt.2018$ci[3,1])
ub_women <- c(out.cent.women.opt.1998$ci[3,2], out.cent.women.opt.2002$ci[3,2]
              , out.cent.women.opt.2006$ci[3,2], out.cent.women.opt.2010$ci[3,2], out.cent.women.opt.2014$ci[3,2]
              , out.cent.women.opt.2018$ci[3,2])

# Calculate Differences 
zOut1998 <- zDif(beta1 = out.cent.men.opt.1998$coef[3,], beta2 = out.cent.women.opt.1998$coef[3,]
                 , se1 = out.cent.men.opt.1998$se[3,], se2 = out.cent.women.opt.1998$se[3,])
zOut2002 <- zDif(beta1 = out.cent.men.opt.2002$coef[3,], beta2 = out.cent.women.opt.2002$coef[3,]
                 , se1 = out.cent.men.opt.2002$se[3,], se2 = out.cent.women.opt.2002$se[3,])
zOut2006 <- zDif(beta1 = out.cent.men.opt.2006$coef[3,], beta2 = out.cent.women.opt.2006$coef[3,]
                 , se1 = out.cent.men.opt.2006$se[3,], se2 = out.cent.women.opt.2006$se[3,])
zOut2010 <- zDif(beta1 = out.cent.men.opt.2010$coef[3,], beta2 = out.cent.women.opt.2010$coef[3,]
                 , se1 = out.cent.men.opt.2010$se[3,], se2 = out.cent.women.opt.2010$se[3,])
zOut2014 <- zDif(beta1 = out.cent.men.opt.2014$coef[3,], beta2 = out.cent.women.opt.2014$coef[3,]
                 , se1 = out.cent.men.opt.2014$se[3,], se2 = out.cent.women.opt.2014$se[3,])
zOut2018 <- zDif(beta1 = out.cent.men.opt.2018$coef[3,], beta2 = out.cent.women.opt.2018$coef[3,]
                 , se1 = out.cent.men.opt.2018$se[3,], se2 = out.cent.women.opt.2018$se[3,])

# FIGURE D.13
# Start the plot:
pdf('./AppendixDResults/figD13.pdf', width = 14)
par(mar = c(5, 5, 2, 2))
plot(x = seq(1998-1, 2018-1, 4), y = betas_men, axes = FALSE, 
     pch = 20, cex = 2, xlim = c(1996, 2020), ylim = c(-3, 3),
     ylab = "RD Effect of Co-Partisan Mayor",
     xlab = "Federal Election", 
     cex.lab = 2)
sapply(1:6, function(i) lines(x = c(seq(1998-1, 2018-1, 4)[i], seq(1998-1, 2018-1, 4)[i]), 
                              y = c(lb_men[i], ub_men[i])))
sapply(1:6, function(i) lines(x = c(seq(1998+1, 2018+1, 4)[i], seq(1998+1, 2018+1, 4)[i]), 
                              y = c(lb_women[i], ub_women[i]), col = 'gray46'))
axis(1, seq(1998, 2018, 4), at = seq(1998, 2018, 4), tick = FALSE, cex.axis = 2)
axis(2, tick = T, cex.axis = 2, las = 2)
segments(x0 = 1995, x1 = 2020.5, y0 = 0, y1 = 0, lty = 3)

# Add differences
# 1998
segments(x0 = 1997, x1 = 1998.5, y0 = betas_men[1], y1 = betas_men[1], lty = 2)
segments(x0 = 1998.5, x1 = 1998.5, y0 = betas_men[1], y1 = betas_women[1], lty = 2)
segments(x0 = 1998.5, x1 = 1999, y0 = betas_women[1], y1 = betas_women[1], lty = 2)
# 2002
segments(x0 = 2001, x1 = 2002.5, y0 = betas_men[2], y1 = betas_men[2], lty = 2)
segments(x0 = 2002.5, x1 = 2002.5, y0 = betas_men[2], y1 = betas_women[2], lty = 2)
segments(x0 = 2002.5, x1 = 2003, y0 = betas_women[2], y1 = betas_women[2], lty = 2)
# 2006
segments(x0 = 2005, x1 = 2006.5, y0 = betas_men[3], y1 = betas_men[3], lty = 2)
segments(x0 = 2006.5, x1 = 2006.5, y0 = betas_men[3], y1 = betas_women[3], lty = 2)
segments(x0 = 2006.5, x1 = 2007, y0 = betas_women[3], y1 = betas_women[3], lty = 2)
# 2010
segments(x0 = 2009, x1 = 2010.5, y0 = betas_men[4], y1 = betas_men[4], lty = 2)
segments(x0 = 2010.5, x1 = 2010.5, y0 = betas_men[4], y1 = betas_women[4], lty = 2)
segments(x0 = 2010.5, x1 = 2011, y0 = betas_women[4], y1 = betas_women[4], lty = 2)
# 2014
segments(x0 = 2013, x1 = 2014.5, y0 = betas_men[5], y1 = betas_men[5], lty = 2)
segments(x0 = 2014.5, x1 = 2014.5, y0 = betas_men[5], y1 = betas_women[5], lty = 2)
segments(x0 = 2014.5, x1 = 2015, y0 = betas_women[5], y1 = betas_women[5], lty = 2)
# 2018
segments(x0 = 2017, x1 = 2018.5, y0 = betas_men[6], y1 = betas_men[6], lty = 2)
segments(x0 = 2018.5, x1 = 2018.5, y0 = betas_men[6], y1 = betas_women[6], lty = 2)
segments(x0 = 2018.5, x1 = 2019, y0 = betas_women[6], y1 = betas_women[6], lty = 2)

# Add points for women's coefficients
points(x = seq(1998+1, 2018+1, 4), y = betas_women, pch = 17, col = 'gray46', cex = 1.5)

# Add information about difference
text(round(zOut1998$difference, 3), y = betas_men[1] + 0.5, x = 1998, cex = 1.5)
text(paste0('(', round(ifelse(zOut1998$zscore > 0, (1 - pnorm(zOut1998$zscore)) * 2, pnorm(zOut1998$zscore) * 2), 3), ')')
     , y = betas_men[1] + 0.2, x = 1998, cex = 1.5)

text(round(zOut2002$difference, 3), y = betas_men[2] + 0.75, x = 2002, cex = 1.5)
text(paste0('(', round(ifelse(zOut2002$zscore > 0, (1 - pnorm(zOut2002$zscore)) * 2, pnorm(zOut2002$zscore) * 2), 3), ')')
     , y = betas_men[2] + 0.45, x = 2002, cex = 1.5)

text(round(zOut2006$difference, 3), y = betas_men[3] + 0.5, x = 2006, cex = 1.5)
text(paste0('(', round(ifelse(zOut2006$zscore > 0, (1 - pnorm(zOut2006$zscore)) * 2, pnorm(zOut2006$zscore) * 2), 3), ')')
     , y = betas_men[3] + 0.2, x = 2006, cex = 1.5)

text(round(zOut2010$difference, 3), y = betas_men[4] + 0.5, x = 2010, cex = 1.5)
text(paste0('(', round(ifelse(zOut2010$zscore > 0, (1 - pnorm(zOut2010$zscore)) * 2, pnorm(zOut2010$zscore) * 2), 3), ')')
     , y = betas_men[4] + 0.2, x = 2010, cex = 1.5)

text(round(zOut2014$difference, 3), y = betas_men[5] + 0.5, x = 2014, cex = 1.5)
text(paste0('(', round(ifelse(zOut2014$zscore > 0, (1 - pnorm(zOut2014$zscore)) * 2, pnorm(zOut2014$zscore) * 2), 3), ')')
     , y = betas_men[5] + 0.2, x = 2014, cex = 1.5)

text(round(zOut2018$difference, 3), y = betas_men[6] + 0.5, x = 2018, cex = 1.5)
text(paste0('(', round(ifelse(zOut2018$zscore > 0, (1 - pnorm(zOut2018$zscore)) * 2, pnorm(zOut2018$zscore) * 2), 3), ')')
     , y = betas_men[6] + 0.2, x = 2018, cex = 1.5)

# Add legend
legend('bottomleft', 
       legend = c('Women', 'Men'), 
       col = c('gray46', 'black'), 
       pch = c(17, 20), 
       lty = c(1, 1),
       bty = 'n',
       cex = 1.5)
dev.off()


###### Manipulation Check
run_varT1 <- d1$run_var[d1$codeyear %in% centParties & d1$ANO_ELEICAO == 1996]
run_varT2 <- d1$run_var[d1$codeyear %in% centParties & d1$ANO_ELEICAO == 2000]
run_varT3 <- d1$run_var[d1$codeyear %in% centParties & d1$ANO_ELEICAO == 2004]
run_varT4 <- d1$run_var[d1$codeyear %in% centParties & d1$ANO_ELEICAO == 2008]
run_varT5 <- d1$run_var[d1$codeyear %in% centParties & d1$ANO_ELEICAO == 2012]
run_varT6 <- d1$run_var[d1$codeyear %in% centParties & d1$ANO_ELEICAO == 2016]
summary(t1 <- rddensity(X = run_varT1, c = 0)) # With massPoints
summary(t2 <- rddensity(X = run_varT2, c = 0)) # With massPoints
summary(t3 <- rddensity(X = run_varT3, c = 0)) # With massPoints
summary(t4 <- rddensity(X = run_varT4, c = 0)) # With massPoints
summary(t5 <- rddensity(X = run_varT5, c = 0)) # With massPoints
summary(t6 <- rddensity(X = run_varT6, c = 0)) # With massPoints
# Result: No sorting

# FIGURES D.14 and D.15
# Density Plot 
pdf('./AppendixDResults/figD14_figD15.pdf', width = 8, height = 8)
par(mar = c(5, 5, 2, 2))
hist(run_varT1[abs(run_varT1)<=(13.039 - 10.873)/2 + 10.873], breaks = 50, # at the bw
     main = "", xlab = "Difference in Vote Share", cex.lab = 1.5, cex.axis = 1.5)
abline(v = 0, col = "red", lwd = 1.5)
text(x = 8.5, y = 22, labels = paste("p-value = ", round(t1$test$p_jk, 4)), cex = 1.25)
par(mar = c(5, 5, 2, 2))
hist(run_varT2[abs(run_varT2)<=(13.039 - 10.873)/2 + 10.873], breaks = 50, # at the bw
     main = "", xlab = "Difference in Vote Share", cex.lab = 1.5, cex.axis = 1.5)
abline(v = 0, col = "red", lwd = 1.5)
text(x = 8.5, y = 30, labels = paste("p-value = ", signif(t2$test$p_jk, 4)), cex = 1.25)
par(mar = c(5, 5, 2, 2))
hist(run_varT3[abs(run_varT3)<=(13.039 - 10.873)/2 + 10.873], breaks = 50, # at the bw
     main = "", xlab = "Difference in Vote Share", cex.lab = 1.5, cex.axis = 1.5)
abline(v = 0, col = "red", lwd = 1.5)
text(x = 8.5, y = 40, labels = paste("p-value = ", signif(t3$test$p_jk, 4)), cex = 1.25)
hist(run_varT4[abs(run_varT4)<=(13.039 - 10.873)/2 + 10.873], breaks = 50, # at the bw
     main = "", xlab = "Difference in Vote Share", cex.lab = 1.5, cex.axis = 1.5)
abline(v = 0, col = "red", lwd = 1.5)
text(x = 8.5, y = 40, labels = paste("p-value = ", round(t4$test$p_jk, 4)), cex = 1.25)
par(mar = c(5, 5, 2, 2))
hist(run_varT5[abs(run_varT5)<=(13.039 - 10.873)/2 + 10.873], breaks = 50, # at the bw
     main = "", xlab = "Difference in Vote Share", cex.lab = 1.5, cex.axis = 1.5)
abline(v = 0, col = "red", lwd = 1.5)
text(x = 8.5, y = 45, labels = paste("p-value = ", signif(t5$test$p_jk, 4)), cex = 1.25)
par(mar = c(5, 5, 2, 2))
hist(run_varT6[abs(run_varT6)<=(13.039 - 10.873)/2 + 10.873], breaks = 50, # at the bw
     main = "", xlab = "Difference in Vote Share", cex.lab = 1.5, cex.axis = 1.5)
abline(v = 0, col = "red", lwd = 1.5)
text(x = 8.5, y = 40, labels = paste("p-value = ", signif(t6$test$p_jk, 4)), cex = 1.25)
dev.off()

#### Lagged Vote Share as the dependent variable
# % of vote share male
d1$perc_party_mun_male_placebo1DF <- d1$perc_party_mun_placebo1DF - d1$perc_party_mun_fem_placebo1DF
# N of male candidates
d1$n_party_Malecandidates_placebo1DF <- d1$n_party_candidates_placebo1DF - d1$n_party_Femcandidates_placebo1DF
# Avg vote share male
d1$avg_placebo1DF <- d1$perc_party_mun_placebo1DF/(d1$n_party_candidates_placebo1DF)
# Avg vote share male
d1$avg_male_placebo1DF <- d1$perc_party_mun_male_placebo1DF/(d1$n_party_Malecandidates_placebo1DF)
# Avg vote share female
d1$avg_fem_placebo1DF <- d1$perc_party_mun_fem_placebo1DF/(d1$n_party_Femcandidates_placebo1DF)

# Find municipalities that should be included (both parties had male and female candidates)
valid <- aggregate(d1$n_party_Malecandidates_placebo1DF > 0 & d1$n_party_Femcandidates_placebo1DF>0 & (!is.na(d1$avg_male_placebo1DF) & !is.na(d1$avg_fem_placebo1DF)), list(d1$codeyear), sum)
valid_place <- na.omit(valid$Group.1[valid$x == 2])
dBalance <- subset(d, codeyear %in% valid_place)

# 2002
out.cent.men.opt.2002 <- rdrobust(y = subset(dBalance,  ANO_ELEICAO.x == 2000)$avg_male_placebo1DF
                                  , x = subset(dBalance, ANO_ELEICAO.x == 2000)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(dBalance,  ANO_ELEICAO.x == 2000)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2002 <- rdrobust(y = subset(dBalance, ANO_ELEICAO.x == 2000)$avg_fem_placebo1DF
                                    , x = subset(dBalance, ANO_ELEICAO.x == 2000)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(dBalance,ANO_ELEICAO.x == 2000)$codeyear
                                    , all = TRUE)

# 2006
out.cent.men.opt.2006 <- rdrobust(y = subset(dBalance,ANO_ELEICAO.x == 2004)$avg_male_placebo1DF
                                  , x = subset(dBalance,ANO_ELEICAO.x == 2004)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(dBalance,ANO_ELEICAO.x == 2004)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2006 <- rdrobust(y = subset(dBalance,ANO_ELEICAO.x == 2004)$avg_fem_placebo1DF
                                    , x = subset(dBalance, ANO_ELEICAO.x == 2004)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(dBalance, ANO_ELEICAO.x == 2004)$codeyear
                                    , all = TRUE)

# 2010
out.cent.men.opt.2010 <- rdrobust(y = subset(dBalance, ANO_ELEICAO.x == 2008)$avg_male_placebo1DF
                                  , x = subset(dBalance, ANO_ELEICAO.x == 2008)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(dBalance, ANO_ELEICAO.x == 2008)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2010 <- rdrobust(y = subset(dBalance, ANO_ELEICAO.x == 2008)$avg_fem_placebo1DF
                                    , x = subset(dBalance, ANO_ELEICAO.x == 2008)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(dBalance, ANO_ELEICAO.x == 2008)$codeyear
                                    , all = TRUE)

# 2014
out.cent.men.opt.2014 <- rdrobust(y = subset(dBalance, ANO_ELEICAO.x == 2012)$avg_male_placebo1DF
                                  , x = subset(dBalance, ANO_ELEICAO.x == 2012)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(dBalance, ANO_ELEICAO.x == 2012)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2014 <- rdrobust(y = subset(dBalance, ANO_ELEICAO.x == 2012)$avg_fem_placebo1DF
                                    , x = subset(dBalance, ANO_ELEICAO.x == 2012)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(dBalance, ANO_ELEICAO.x == 2012)$codeyear
                                    , all = TRUE)

# 2018
out.cent.men.opt.2018 <- rdrobust(y = subset(dBalance, ANO_ELEICAO.x == 2016)$avg_male_placebo1DF
                                  , x = subset(dBalance, ANO_ELEICAO.x == 2016)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(dBalance, ANO_ELEICAO.x == 2016)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2018 <- rdrobust(y = subset(dBalance, ANO_ELEICAO.x == 2016)$avg_fem_placebo1DF
                                    , x = subset(dBalance, ANO_ELEICAO.x == 2016)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(dBalance, ANO_ELEICAO.x == 2016)$codeyear
                                    , all = TRUE)

# Get point estimates, se, and ci for men
betas_men <- c(out.cent.men.opt.2002$coef[3,]
               , out.cent.men.opt.2006$coef[3,], out.cent.men.opt.2010$coef[3,], out.cent.men.opt.2014$coef[3,]
               , out.cent.men.opt.2018$coef[3,])
se_men <- c(out.cent.men.opt.2002$se[3,]
            , out.cent.men.opt.2006$se[3,], out.cent.men.opt.2010$se[3,], out.cent.men.opt.2014$se[3,]
            , out.cent.men.opt.2018$se[3,])
lb_men <- c(out.cent.men.opt.2002$ci[3,1]
            , out.cent.men.opt.2006$ci[3,1], out.cent.men.opt.2010$ci[3,1], out.cent.men.opt.2014$ci[3,1]
            , out.cent.men.opt.2018$ci[3,1])
ub_men <- c(out.cent.men.opt.2002$ci[3,2]
            , out.cent.men.opt.2006$ci[3,2], out.cent.men.opt.2010$ci[3,2], out.cent.men.opt.2014$ci[3,2]
            , out.cent.men.opt.2018$ci[3,2])

# Get point estimates, se, and ci for women
betas_women <- c(out.cent.women.opt.2002$coef[3,], out.cent.women.opt.2006$coef[3,]
                 , out.cent.women.opt.2010$coef[3,], out.cent.women.opt.2014$coef[3,]
                 , out.cent.women.opt.2018$coef[3,])
se_women <- c(out.cent.women.opt.2002$se[3,]
              , out.cent.women.opt.2006$se[3,], out.cent.women.opt.2010$se[3,], out.cent.women.opt.2014$se[3,]
              , out.cent.women.opt.2018$se[3,])
lb_women <- c(out.cent.women.opt.2002$ci[3,1]
              , out.cent.women.opt.2006$ci[3,1], out.cent.women.opt.2010$ci[3,1], out.cent.women.opt.2014$ci[3,1]
              , out.cent.women.opt.2018$ci[3,1])
ub_women <- c(out.cent.women.opt.2002$ci[3,2]
              , out.cent.women.opt.2006$ci[3,2], out.cent.women.opt.2010$ci[3,2], out.cent.women.opt.2014$ci[3,2]
              , out.cent.women.opt.2018$ci[3,2])

# Calculate Differences 
zOut <- list(
  zDif(beta1 = out.cent.men.opt.2002$coef[3,], beta2 = out.cent.women.opt.2002$coef[3,]
       , se1 = out.cent.men.opt.2002$se[3,], se2 = out.cent.women.opt.2002$se[3,]),
  zDif(beta1 = out.cent.men.opt.2006$coef[3,], beta2 = out.cent.women.opt.2006$coef[3,]
       , se1 = out.cent.men.opt.2006$se[3,], se2 = out.cent.women.opt.2006$se[3,]),
  zDif(beta1 = out.cent.men.opt.2010$coef[3,], beta2 = out.cent.women.opt.2010$coef[3,]
       , se1 = out.cent.men.opt.2010$se[3,], se2 = out.cent.women.opt.2010$se[3,]),
  zDif(beta1 = out.cent.men.opt.2014$coef[3,], beta2 = out.cent.women.opt.2014$coef[3,]
       , se1 = out.cent.men.opt.2014$se[3,], se2 = out.cent.women.opt.2014$se[3,]),
  zDif(beta1 = out.cent.men.opt.2018$coef[3,], beta2 = out.cent.women.opt.2018$coef[3,]
       , se1 = out.cent.men.opt.2018$se[3,], se2 = out.cent.women.opt.2018$se[3,])
)


# FIGURE D.16
pdf('./AppendixDResults/figD16.pdf', width = 14)
byYearH1(betas_men = betas_men, betas_women = betas_women, lb_men = lb_men, ub_men = ub_men, 
         lb_women = lb_women, ub_women = ub_women, zout = zOut, ylims = c(-5, 5), betaPos = 0.30)
dev.off()


##### Sensitivity
# FIGURE D.17
genSensitivityH1(yearElection = 1996, yearPlot = 1998, file="./AppendixDResults/figD17.pdf")

# FIGURE D.18
genSensitivityH1(yearElection = 2000, yearPlot = 2002, file="./AppendixDResults/figD18.pdf")

# FIGURE D.19
genSensitivityH1(yearElection = 2004, yearPlot = 2006, file="./AppendixDResults/figD19.pdf")

# FIGURE D.20
genSensitivityH1(yearElection = 2008, yearPlot = 2010, file="./AppendixDResults/figD20.pdf")

# FIGURE D.21
genSensitivityH1(yearElection = 2012, yearPlot = 2014, file="./AppendixDResults/figD21.pdf")

# FIGURE D.22
genSensitivityH1(yearElection = 2016, yearPlot = 2018, file="./AppendixDResults/figD22.pdf")

# TABLES D.14 and D.15
##### Different Polynomials

sink("./AppendixDResults/tableD14.txt")
print("TABLE D.14: RD Effect of Co-Partisan Mayor on Average Vote Share (Different Polynomials)—
Brazilian Federal Elections, 1998-2006")

polyH1(yearElection = 1996, p = 2)
polyH1(yearElection = 1996, p = 3)
polyH1(yearElection = 1996, p = 4)

polyH1(yearElection = 2000, p = 2)
polyH1(yearElection = 2000, p = 3)
polyH1(yearElection = 2000, p = 4)

polyH1(yearElection = 2004, p = 2)
polyH1(yearElection = 2004, p = 3)
polyH1(yearElection = 2004, p = 4)
sink()

sink("./AppendixDResults/tableD15.txt")
print("TABLE D.15: RD Effect of Co-Partisan Mayor on Average Vote Share (Different Polynomials)—
Brazilian Federal Elections, 2010-2018")

polyH1(yearElection = 2008, p = 2)
polyH1(yearElection = 2008, p = 3)
polyH1(yearElection = 2008, p = 4)

polyH1(yearElection = 2012, p = 2)
polyH1(yearElection = 2012, p = 3)
polyH1(yearElection = 2012, p = 4)

polyH1(yearElection = 2016, p = 2)
polyH1(yearElection = 2016, p = 3)
polyH1(yearElection = 2016, p = 4)
sink()


##### Previously Ran

# N of male candidates
d1$n_party_candidates_male_candBefore_DF <- d1$n_party_candidates_candBefore_DF - d1$n_party_Femcandidates_candBefore_DF
valid <- aggregate(d1$n_party_candidates_male_candBefore_DF > 0 & d1$n_party_Femcandidates_candBefore_DF>0, list(d1$codeyear), sum)
valid_ranB <- na.omit(valid$Group.1[valid$x == 2])
# % of vote share male
d1$perc_party_mun_male_candBeforeDF <- d1$perc_party_mun_candBefore_DF - d1$perc_party_mun_fem_candBefore_DF
# Avg vote share all
d1$avg_candBefore_DF <- d1$perc_party_mun_candBefore_DF/d1$n_party_candidates_candBefore_DF
# Avg vote share male
d1$avg_male_candBefore_DF <- d1$perc_party_mun_male_candBeforeDF/(d1$n_party_candidates_candBefore_DF - d1$n_party_Femcandidates_candBefore_DF)
# Avg vote share female
d1$avg_fem_candBefore_DF <- d1$perc_party_mun_fem_candBefore_DF/(d1$n_party_Femcandidates_candBefore_DF)
# Subset data
d2 <- subset(d1, codeyear %in% valid_ranB)

# 2002
out.cent.men.opt.2002 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2000)$avg_male_candBefore_DF
                                  , x = subset(d2, ANO_ELEICAO == 2000)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d2, ANO_ELEICAO == 2000)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2002 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2000)$avg_fem_candBefore_DF
                                    , x = subset(d2, ANO_ELEICAO == 2000)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d2, ANO_ELEICAO == 2000)$codeyear
                                    , all = TRUE)

# 2006
out.cent.men.opt.2006 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2004)$avg_male_candBefore_DF
                                  , x = subset(d2, ANO_ELEICAO == 2004)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d2, ANO_ELEICAO == 2004)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2006 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2004)$avg_fem_candBefore_DF
                                    , x = subset(d2, ANO_ELEICAO == 2004)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d2, ANO_ELEICAO == 2004)$codeyear
                                    , all = TRUE)

# 2010
out.cent.men.opt.2010 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2008)$avg_male_candBefore_DF
                                  , x = subset(d2, ANO_ELEICAO == 2008)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d2, ANO_ELEICAO == 2008)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2010 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2008)$avg_fem_candBefore_DF
                                    , x = subset(d2, ANO_ELEICAO == 2008)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d2, ANO_ELEICAO == 2008)$codeyear
                                    , all = TRUE)

# 2014
out.cent.men.opt.2014 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2012)$avg_male_candBefore_DF
                                  , x = subset(d2, ANO_ELEICAO == 2012)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d2, ANO_ELEICAO == 2012)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2014 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2012)$avg_fem_candBefore_DF
                                    , x = subset(d2, ANO_ELEICAO == 2012)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d2, ANO_ELEICAO == 2012)$codeyear
                                    , all = TRUE)

# 2018
out.cent.men.opt.2018 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2016)$avg_male_candBefore_DF
                                  , x = subset(d2, ANO_ELEICAO == 2016)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d2, ANO_ELEICAO == 2016)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2018 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2016)$avg_fem_candBefore_DF
                                    , x = subset(d2, ANO_ELEICAO == 2016)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d2, ANO_ELEICAO == 2016)$codeyear
                                    , all = TRUE)

# Get point estimates, se, and ci for men
betas_men <- c(out.cent.men.opt.2002$coef[3,]
               , out.cent.men.opt.2006$coef[3,], out.cent.men.opt.2010$coef[3,], out.cent.men.opt.2014$coef[3,]
               , out.cent.men.opt.2018$coef[3,])
se_men <- c(out.cent.men.opt.2002$se[3,]
            , out.cent.men.opt.2006$se[3,], out.cent.men.opt.2010$se[3,], out.cent.men.opt.2014$se[3,]
            , out.cent.men.opt.2018$se[3,])
lb_men <- c(out.cent.men.opt.2002$ci[3,1]
            , out.cent.men.opt.2006$ci[3,1], out.cent.men.opt.2010$ci[3,1], out.cent.men.opt.2014$ci[3,1]
            , out.cent.men.opt.2018$ci[3,1])
ub_men <- c(out.cent.men.opt.2002$ci[3,2]
            , out.cent.men.opt.2006$ci[3,2], out.cent.men.opt.2010$ci[3,2], out.cent.men.opt.2014$ci[3,2]
            , out.cent.men.opt.2018$ci[3,2])

# Get point estimates, se, and ci for women
betas_women <- c(out.cent.women.opt.2002$coef[3,], out.cent.women.opt.2006$coef[3,]
                 , out.cent.women.opt.2010$coef[3,], out.cent.women.opt.2014$coef[3,]
                 , out.cent.women.opt.2018$coef[3,])
se_women <- c(out.cent.women.opt.2002$se[3,]
              , out.cent.women.opt.2006$se[3,], out.cent.women.opt.2010$se[3,], out.cent.women.opt.2014$se[3,]
              , out.cent.women.opt.2018$se[3,])
lb_women <- c(out.cent.women.opt.2002$ci[3,1]
              , out.cent.women.opt.2006$ci[3,1], out.cent.women.opt.2010$ci[3,1], out.cent.women.opt.2014$ci[3,1]
              , out.cent.women.opt.2018$ci[3,1])
ub_women <- c(out.cent.women.opt.2002$ci[3,2]
              , out.cent.women.opt.2006$ci[3,2], out.cent.women.opt.2010$ci[3,2], out.cent.women.opt.2014$ci[3,2]
              , out.cent.women.opt.2018$ci[3,2])

# Calculate Differences 
zOut <- list(
  zDif(beta1 = out.cent.men.opt.2002$coef[3,], beta2 = out.cent.women.opt.2002$coef[3,]
       , se1 = out.cent.men.opt.2002$se[3,], se2 = out.cent.women.opt.2002$se[3,]),
  zDif(beta1 = out.cent.men.opt.2006$coef[3,], beta2 = out.cent.women.opt.2006$coef[3,]
       , se1 = out.cent.men.opt.2006$se[3,], se2 = out.cent.women.opt.2006$se[3,]),
  zDif(beta1 = out.cent.men.opt.2010$coef[3,], beta2 = out.cent.women.opt.2010$coef[3,]
       , se1 = out.cent.men.opt.2010$se[3,], se2 = out.cent.women.opt.2010$se[3,]),
  zDif(beta1 = out.cent.men.opt.2014$coef[3,], beta2 = out.cent.women.opt.2014$coef[3,]
       , se1 = out.cent.men.opt.2014$se[3,], se2 = out.cent.women.opt.2014$se[3,]),
  zDif(beta1 = out.cent.men.opt.2018$coef[3,], beta2 = out.cent.women.opt.2018$coef[3,]
       , se1 = out.cent.men.opt.2018$se[3,], se2 = out.cent.women.opt.2018$se[3,])
)

# FIGURE D.23
pdf('./AppendixDResults/figD23.pdf', width = 14)
byYearH1(betas_men = betas_men, betas_women = betas_women, lb_men = lb_men, ub_men = ub_men, lb_women = lb_women, ub_women = ub_women, zout = zOut)
dev.off()


##### Previously Won

# Women Candidate who won before
d1$perc_party_mun_fem_wonBefore_DF <- (d1$total_party_mun_fem_wonBefore_DF/d1$total_mun_wonBefore_DF) * 100
# # % of vote share male
d1$perc_party_mun_male_wonBefore_DF <- d1$perc_party_mun_wonBefore_DF - d1$perc_party_mun_fem_wonBefore_DF
# # N of male candidates
d1$n_party_candidates_male_wonBefore_DF <- d1$n_party_candidates_wonBefore_DF - d1$n_party_Femcandidates_wonBefore_DF
# # Avg vote share 
d1$avg_wonBefore_DF <- d1$perc_party_mun_wonBefore_DF/d1$n_party_candidates_wonBefore_DF
# # Avg vote share male
d1$avg_male_wonBefore_DF <- d1$perc_party_mun_male_wonBefore_DF/d1$n_party_candidates_male_wonBefore_DF
# # Avg vote share female
d1$avg_fem_wonBefore_DF <- d1$perc_party_mun_fem_wonBefore_DF/(d1$n_party_Femcandidates_wonBefore_DF)

# # Find municipalities that should be included
valid <- aggregate(d1$n_party_candidates_male_wonBefore_DF > 0 & d1$n_party_Femcandidates_wonBefore_DF>0, list(d1$codeyear), sum)
valid_wonB <- na.omit(valid$Group.1[valid$x == 2])

# Subset data
d2 <- subset(d1, codeyear %in% valid_wonB)

# 2002
out.cent.men.opt.2002 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2000)$avg_male_wonBefore_DF
                                  , x = subset(d2, ANO_ELEICAO == 2000)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d2, ANO_ELEICAO == 2000)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2002 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2000)$avg_fem_wonBefore_DF
                                    , x = subset(d2, ANO_ELEICAO == 2000)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d2, ANO_ELEICAO == 2000)$codeyear
                                    , all = TRUE)

# 2006
out.cent.men.opt.2006 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2004)$avg_male_wonBefore_DF
                                  , x = subset(d2, ANO_ELEICAO == 2004)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d2, ANO_ELEICAO == 2004)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2006 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2004)$avg_fem_wonBefore_DF
                                    , x = subset(d2, ANO_ELEICAO == 2004)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d2, ANO_ELEICAO == 2004)$codeyear
                                    , all = TRUE)

# 2010
out.cent.men.opt.2010 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2008)$avg_male_wonBefore_DF
                                  , x = subset(d2, ANO_ELEICAO == 2008)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d2, ANO_ELEICAO == 2008)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2010 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2008)$avg_fem_wonBefore_DF
                                    , x = subset(d2, ANO_ELEICAO == 2008)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d2, ANO_ELEICAO == 2008)$codeyear
                                    , all = TRUE)

# 2014
out.cent.men.opt.2014 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2012)$avg_male_wonBefore_DF
                                  , x = subset(d2, ANO_ELEICAO == 2012)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d2, ANO_ELEICAO == 2012)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2014 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2012)$avg_fem_wonBefore_DF
                                    , x = subset(d2, ANO_ELEICAO == 2012)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d2, ANO_ELEICAO == 2012)$codeyear
                                    , all = TRUE)

# 2018
out.cent.men.opt.2018 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2016)$avg_male_wonBefore_DF
                                  , x = subset(d2, ANO_ELEICAO == 2016)$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d2, ANO_ELEICAO == 2016)$codeyear
                                  , all = TRUE)
out.cent.women.opt.2018 <- rdrobust(y = subset(d2, ANO_ELEICAO == 2016)$avg_fem_wonBefore_DF
                                    , x = subset(d2, ANO_ELEICAO == 2016)$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d2, ANO_ELEICAO == 2016)$codeyear
                                    , all = TRUE)

# Get point estimates, se, and ci for men
betas_men <- c(out.cent.men.opt.2002$coef[3,]
               , out.cent.men.opt.2006$coef[3,], out.cent.men.opt.2010$coef[3,], out.cent.men.opt.2014$coef[3,]
               , out.cent.men.opt.2018$coef[3,])
se_men <- c(out.cent.men.opt.2002$se[3,]
            , out.cent.men.opt.2006$se[3,], out.cent.men.opt.2010$se[3,], out.cent.men.opt.2014$se[3,]
            , out.cent.men.opt.2018$se[3,])
lb_men <- c(out.cent.men.opt.2002$ci[3,1]
            , out.cent.men.opt.2006$ci[3,1], out.cent.men.opt.2010$ci[3,1], out.cent.men.opt.2014$ci[3,1]
            , out.cent.men.opt.2018$ci[3,1])
ub_men <- c(out.cent.men.opt.2002$ci[3,2]
            , out.cent.men.opt.2006$ci[3,2], out.cent.men.opt.2010$ci[3,2], out.cent.men.opt.2014$ci[3,2]
            , out.cent.men.opt.2018$ci[3,2])

# Get point estimates, se, and ci for women
betas_women <- c(out.cent.women.opt.2002$coef[3,], out.cent.women.opt.2006$coef[3,]
                 , out.cent.women.opt.2010$coef[3,], out.cent.women.opt.2014$coef[3,]
                 , out.cent.women.opt.2018$coef[3,])
se_women <- c(out.cent.women.opt.2002$se[3,]
              , out.cent.women.opt.2006$se[3,], out.cent.women.opt.2010$se[3,], out.cent.women.opt.2014$se[3,]
              , out.cent.women.opt.2018$se[3,])
lb_women <- c(out.cent.women.opt.2002$ci[3,1]
              , out.cent.women.opt.2006$ci[3,1], out.cent.women.opt.2010$ci[3,1], out.cent.women.opt.2014$ci[3,1]
              , out.cent.women.opt.2018$ci[3,1])
ub_women <- c(out.cent.women.opt.2002$ci[3,2]
              , out.cent.women.opt.2006$ci[3,2], out.cent.women.opt.2010$ci[3,2], out.cent.women.opt.2014$ci[3,2]
              , out.cent.women.opt.2018$ci[3,2])

# Calculate Differences 
zOut <- list(
  zDif(beta1 = out.cent.men.opt.2002$coef[3,], beta2 = out.cent.women.opt.2002$coef[3,]
       , se1 = out.cent.men.opt.2002$se[3,], se2 = out.cent.women.opt.2002$se[3,]),
  zDif(beta1 = out.cent.men.opt.2006$coef[3,], beta2 = out.cent.women.opt.2006$coef[3,]
       , se1 = out.cent.men.opt.2006$se[3,], se2 = out.cent.women.opt.2006$se[3,]),
  zDif(beta1 = out.cent.men.opt.2010$coef[3,], beta2 = out.cent.women.opt.2010$coef[3,]
       , se1 = out.cent.men.opt.2010$se[3,], se2 = out.cent.women.opt.2010$se[3,]),
  zDif(beta1 = out.cent.men.opt.2014$coef[3,], beta2 = out.cent.women.opt.2014$coef[3,]
       , se1 = out.cent.men.opt.2014$se[3,], se2 = out.cent.women.opt.2014$se[3,]),
  zDif(beta1 = out.cent.men.opt.2018$coef[3,], beta2 = out.cent.women.opt.2018$coef[3,]
       , se1 = out.cent.men.opt.2018$se[3,], se2 = out.cent.women.opt.2018$se[3,])
)


# FIGURE D.24
pdf('./AppendixDResults/figD24.pdf', width = 14)
byYearH1(betas_men = betas_men, betas_women = betas_women, lb_men = lb_men, ub_men = ub_men, 
         lb_women = lb_women, ub_women = ub_women, zout = zOut, ylims = c(-5, 5), betaPos = 0.25)
dev.off()

########################################
####### RD using Linear Estimator ######
########################################
# Create Dataset
dmale <- subset(d1, select = c('avg_male_DF', 'position', 'run_var', 'codeyear', 'mayorFem', 'ANO_ELEICAO'))
dmale$men <- 1
dfem <- subset(d1, select = c('avg_fem_DF', 'position', 'run_var', 'codeyear', 'mayorFem', 'ANO_ELEICAO'))
dfem$men <- 0
names(dmale)[1] <- 'avg_DF'
names(dfem)[1] <- 'avg_DF'
dall <- rbind(dmale, dfem)
dall$winner <- as.numeric(dall$run_var > 0)
d2 <- subset(dall, codeyear %in% centParties)
# Run Models
out1 <- linerEstimator(1996)
out2 <- linerEstimator(2000)
out3 <- linerEstimator(2004)
out4 <- linerEstimator(2008)
out5 <- linerEstimator(2012)
out6 <- linerEstimator(2016)

# Generate plots
# FIGURE D.25
pdf('./AppendixDResults/figD25.pdf')
plotTtest(betas = out1[[1]], lb = out1[[2]], ub = out1[[3]], dvName = 'Federal Deputy - Men', yub = 4, ylb = -4)
plotTtest(betas = out1[[4]], lb = out1[[5]], ub = out1[[6]], dvName = 'Federal Deputy - Women', yub = 4, ylb = -4)
plotTtest(betas = out1[[7]], lb = out1[[8]], ub = out1[[9]], dvName = 'Federal Deputy - Difference', yub = 4, ylb = -4)
dev.off()
# FIGURE D.26
pdf('./AppendixDResults/figD26.pdf')
plotTtest(betas = out2[[1]], lb = out2[[2]], ub = out2[[3]], dvName = 'Federal Deputy - Men', yub = 4, ylb = -4)
plotTtest(betas = out2[[4]], lb = out2[[5]], ub = out2[[6]], dvName = 'Federal Deputy - Women', yub = 4, ylb = -4)
plotTtest(betas = out2[[7]], lb = out2[[8]], ub = out2[[9]], dvName = 'Federal Deputy - Difference', yub = 4, ylb = -4)
dev.off()
# FIGURE D.27
pdf('./AppendixDResults/figD27.pdf')
plotTtest(betas = out3[[1]], lb = out3[[2]], ub = out3[[3]], dvName = 'Federal Deputy - Men', yub = 4, ylb = -4)
plotTtest(betas = out3[[4]], lb = out3[[5]], ub = out3[[6]], dvName = 'Federal Deputy - Women', yub = 4, ylb = -4)
plotTtest(betas = out3[[7]], lb = out3[[8]], ub = out3[[9]], dvName = 'Federal Deputy - Difference', yub = 4, ylb = -4)
dev.off()
# FIGURE D.28
pdf('./AppendixDResults/figD28.pdf')
plotTtest(betas = out4[[1]], lb = out4[[2]], ub = out4[[3]], dvName = 'Federal Deputy - Men', yub = 4, ylb = -4)
plotTtest(betas = out4[[4]], lb = out4[[5]], ub = out4[[6]], dvName = 'Federal Deputy - Women', yub = 4, ylb = -4)
plotTtest(betas = out4[[7]], lb = out4[[8]], ub = out4[[9]], dvName = 'Federal Deputy - Difference', yub = 4, ylb = -4)
dev.off()
# FIGURE D.29
pdf('./AppendixDResults/figD29.pdf')
plotTtest(betas = out5[[1]], lb = out5[[2]], ub = out5[[3]], dvName = 'Federal Deputy - Men', yub = 4, ylb = -4)
plotTtest(betas = out5[[4]], lb = out5[[5]], ub = out5[[6]], dvName = 'Federal Deputy - Women', yub = 4, ylb = -4)
plotTtest(betas = out5[[7]], lb = out5[[8]], ub = out5[[9]], dvName = 'Federal Deputy - Difference', yub = 4, ylb = -4)
dev.off()
# FIGURE D.30
pdf('./AppendixDResults/figD30.pdf')
plotTtest(betas = out6[[1]], lb = out6[[2]], ub = out6[[3]], dvName = 'Federal Deputy - Men', yub = 4, ylb = -4)
plotTtest(betas = out6[[4]], lb = out6[[5]], ub = out6[[6]], dvName = 'Federal Deputy - Women', yub = 4, ylb = -4)
plotTtest(betas = out6[[7]], lb = out6[[8]], ub = out6[[9]], dvName = 'Federal Deputy - Difference', yub = 4, ylb = -4)
dev.off()

##############################################################################
################# Appendix D.10: Subsample Analysis by Party #################
##############################################################################

# FIGURES D.31, D.32, D.33, D.34, D.35, and D.36

# We are unable to run the model for PC do B due to the limited sample size; therefore, the PCdoB is the only centralized party
# for which we do not produce a subsample analysis

# PT (centralized): FIGURE D.31
figure3ByParty(partyName = 'PT', file = "./AppendixDResults/figD31.pdf")

# PSDB (centralized): FIGURE D.32
figure3ByParty(partyName = 'PSDB', file = "./AppendixDResults/figD32.pdf")

# PPS (centralized): FIGURE D.33
figure3ByParty(partyName = 'PPS', file = "./AppendixDResults/figD33.pdf")

# DEM/PFL (centralized): FIGURE D.34
figure3ByParty(partyName = c('PFL', 'DEM'), file = "./AppendixDResults/figD34.pdf")

# PL/PR (centralized): FIGURE D.35
figure3ByParty(partyName = c('PL', 'PR'), file = "./AppendixDResults/figD35.pdf")

# PDT (centralized): FIGURE D.36
figure3ByParty(partyName = 'PDT', file = "./AppendixDResults/figD36.pdf")


##############################################################################
####################### Appendix D.11: Decentralized Parties #################
##############################################################################

# Weak
out.cent.men.opt.weak <- rdrobust(y = subset(d_d1,  period == 'weak')$avg_male_DF
                                  , x = subset(d_d1,  period == 'weak')$run_var
                                  , p = 1
                                  , c = 0	
                                  , cluster = subset(d_d1, period == 'weak')$codeyear
                                  , all = TRUE)
out.cent.women.opt.weak <- rdrobust(y = subset(d_d1, period == 'weak')$avg_fem_DF
                                    , x = subset(d_d1,  period == 'weak')$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d_d1, period == 'weak')$codeyear
                                    , all = TRUE)

# Strong
out.cent.men.opt.strong <- rdrobust(y = subset(d_d1, period == 'strong')$avg_male_DF
                                    , x = subset(d_d1, period == 'strong')$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d_d1, period == 'strong')$codeyear
                                    , all = TRUE)
out.cent.women.opt.strong <- rdrobust(y = subset(d_d1, period == 'strong')$avg_fem_DF
                                      , x = subset(d_d1, period == 'strong')$run_var
                                      , p = 1
                                      , c = 0	
                                      , cluster = subset(d_d1, period == 'strong')$codeyear
                                      , all = TRUE)

# Moderate
out.cent.men.opt.moderate <- rdrobust(y = subset(d_d1, period == 'moderate')$avg_male_DF
                                      , x = subset(d_d1, period == 'moderate')$run_var
                                      , p = 1
                                      , c = 0	
                                      , cluster = subset(d_d1, period == 'moderate')$codeyear
                                      , all = TRUE)
out.cent.women.opt.moderate <- rdrobust(y = subset(d_d1, period == 'moderate')$avg_fem_DF
                                        , x = subset(d_d1, period == 'moderate')$run_var
                                        , p = 1
                                        , c = 0	
                                        , cluster = subset(d_d1, period == 'moderate')$codeyear
                                        , all = TRUE)

# Get point estimates, se, and ci for men
betas_men <- c(out.cent.men.opt.weak$coef[3,], out.cent.men.opt.strong$coef[3,]
               , out.cent.men.opt.moderate$coef[3,])
se_men <- c(out.cent.men.opt.weak$se[3,], out.cent.men.opt.strong$se[3,]
            , out.cent.men.opt.moderate$se[3,])
lb_men <- c(out.cent.men.opt.weak$ci[3,1], out.cent.men.opt.strong$ci[3,1]
            , out.cent.men.opt.moderate$ci[3,1])
ub_men <- c(out.cent.men.opt.weak$ci[3,2], out.cent.men.opt.strong$ci[3,2]
            , out.cent.men.opt.moderate$ci[3,2])

# Get point estimates, se, and ci for women
betas_women <- c(out.cent.women.opt.weak$coef[3,], out.cent.women.opt.strong$coef[3,]
                 , out.cent.women.opt.moderate$coef[3,])
se_women <- c(out.cent.women.opt.weak$se[3,], out.cent.women.opt.strong$se[3,]
              , out.cent.women.opt.moderate$se[3,])
lb_women <- c(out.cent.women.opt.weak$ci[3,1], out.cent.women.opt.strong$ci[3,1]
              , out.cent.women.opt.moderate$ci[3,1])
ub_women <- c(out.cent.women.opt.weak$ci[3,2], out.cent.women.opt.strong$ci[3,2]
              , out.cent.women.opt.moderate$ci[3,2])

# Calculate Differences 
zOutWeak <- zDif(beta1 = out.cent.men.opt.weak$coef[3,], beta2 = out.cent.women.opt.weak$coef[3,]
                 , se1 = out.cent.men.opt.weak$se[3,], se2 = out.cent.women.opt.weak$se[3,])
zOutStrong <- zDif(beta1 = out.cent.men.opt.strong$coef[3,], beta2 = out.cent.women.opt.strong$coef[3,]
                   , se1 = out.cent.men.opt.strong$se[3,], se2 = out.cent.women.opt.strong$se[3,])
zOutModerate <- zDif(beta1 = out.cent.men.opt.moderate$coef[3,], beta2 = out.cent.women.opt.moderate$coef[3,]
                     , se1 = out.cent.men.opt.moderate$se[3,], se2 = out.cent.women.opt.moderate$se[3,])


# Start the plot:
pdf('./AppendixDResults/figD37.pdf', width = 14)
par(mar = c(5, 5, 2, 2))
plot(x = seq(1-0.1, 3-0.1, 1), y = betas_men, axes = FALSE, 
     pch = 20, cex = 2, xlim = c(0.75, 3.25), ylim = c(-3, 3),
     ylab = "RD Effect of Co-Partisan Mayor",
     xlab = "Electoral Rule Regime", 
     cex.lab = 2)
sapply(1:3, function(i) lines(x = c(seq(1-0.1, 3-0.1, 1)[i], seq(1-0.1, 3-0.1, 1)[i]), 
                              y = c(lb_men[i], ub_men[i])))
sapply(1:3, function(i) lines(x = c(seq(1+0.1, 3+0.1, 1)[i], seq(1+0.1, 3+0.1, 1)[i]), 
                              y = c(lb_women[i], ub_women[i]), col = 'gray46'))
axis(1, at = seq(1, 3, 1), label = c('Weak\n(1998-2006)', 'Strong\n(2010-2014)', 'Moderate\n(2018)'), tick = FALSE, cex.axis = 2)
axis(2, tick = T, cex.axis = 2, las = 2)
segments(x0 = 0, x1 = 3.5, y0 = 0, y1 = 0, lty = 3)

# Add differences
# Weak
segments(x0 = 0.9, x1 = 1.05, y0 = betas_men[1], y1 = betas_men[1], lty = 2)
segments(x0 = 1.05, x1 = 1.05, y0 = betas_men[1], y1 = betas_women[1], lty = 2)
segments(x0 = 1.05, x1 = 1.1, y0 = betas_women[1], y1 = betas_women[1], lty = 2)
# Strong
segments(x0 = 1.9, x1 = 2.05, y0 = betas_men[2], y1 = betas_men[2], lty = 2)
segments(x0 = 2.05, x1 = 2.05, y0 = betas_men[2], y1 = betas_women[2], lty = 2)
segments(x0 = 2.05, x1 = 2.1, y0 = betas_women[2], y1 = betas_women[2], lty = 2)
# Moderate
segments(x0 = 2.9, x1 = 3.05, y0 = betas_men[3], y1 = betas_men[3], lty = 2)
segments(x0 = 3.05, x1 = 3.05, y0 = betas_men[3], y1 = betas_women[3], lty = 2)
segments(x0 = 3.05, x1 = 3.1, y0 = betas_women[3], y1 = betas_women[3], lty = 2)

# Add points for women's coefficients
points(x = seq(1+0.1, 3+0.1, 1), y = betas_women, pch = 17, col = 'gray46', cex = 1.5)

# Add information about difference
text(round(zOutWeak$difference, 3), y = betas_men[1] + 0.5, x = 1, cex = 1.5)
text(paste0('(', round(ifelse(zOutWeak$zscore > 0, (1 - pnorm(zOutWeak$zscore)) * 2, pnorm(zOutWeak$zscore) * 2), 3), ')')
     , y = betas_men[1] + 0.2, x = 1, cex = 1.5)

text(round(zOutStrong$difference, 3), y = betas_men[2] + 0.75, x = 2, cex = 1.5)
text(paste0('(', ifelse(round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3) == 0, "0.000", round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3)), ')')
     , y = betas_men[2] + 0.45, x = 2, cex = 1.5)

text(round(zOutModerate$difference, 3), y = betas_men[3] + 0.5, x = 3, cex = 1.5)
text(paste0('(', round(ifelse(zOutModerate$zscore > 0, (1 - pnorm(zOutModerate$zscore)) * 2, pnorm(zOutModerate$zscore) * 2), 3), ')')
     , y = betas_men[3] + 0.2, x = 3, cex = 1.5)

# Add legend
legend(x = 0.7, y = -1.5, 
       legend = c('Women', 'Men'), 
       col = c('gray46', 'black'), 
       pch = c(17, 20), 
       lty = c(1, 1),
       bty = 'n',
       cex = 1.5)
dev.off()
